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In this study, we analyze the statistics of both individual inertial particles and inertial 
particle pairs in direct numerical simulations of homogeneous isotropic turbulence in 
the absence of gravity. The effect of the Taylor microscale Reynolds number R\ on the 
particle statistics is examined over the largest range to date (from R\ = 88 — 597), at 
small, intermediate, and large Kolmogorov-scale Stokes numbers St. We first explore the 
effect of preferential sampling on the single-particle statistics and find that low-S't inertial 
particles are ejected from both vortex tubes and vortex sheets (the latter becoming 
increasingly prevalent at higher Reynolds numbers) and preferentially accumulate in 
regions of irrotational dissipation. We use this understanding of preferential sampling to 
provide a physical explanation for many of the trends in the particle velocity gradients, 
kinetic energies, and accelerations at low St, which are well-represented by the model 
of Chun et al. (2005). As St increases, inertial filtering effects become more important, 
causing the particle kinetic energies and accelerations to decrease. The effect of inertial 
filtering on the particle kinetic energies and accelerations diminishes with increasing 
Reynolds number and is well-captured by the models of Abrahamson (1975) and Zaichik 
& Alipchenkov (2008), respectively. 

We then consider particle-pair statistics, and focus our attention on the relative ve¬ 
locities and radial distribution functions (RDFs) of the particles, with the aim of un¬ 
derstanding the underlying physical mechanisms contributing to particle collisions. The 
relative velocity statistics indicate that preferential-sampling effects are important for 
St < 0.1 and that path-history/non-local effects become increasingly important for 
St > 0.2. While higher-order relative velocity statistics are influenced by the increased 
intermittency of the turbulence at high Reynolds numbers, the lower-order relative ve¬ 
locity statistics are only weakly sensitive to changes in Reynolds number at low St. The 
Reynolds-number trends in these quantities at intermediate and large St are explained 
based on the influence of the available flow scales on the path-history and inertial filtering 
effects. We find that the RDFs peak near St of order unity, that they exhibit power-law 
scaling for low and intermediate St, and that they are largely independent of Reynolds 
number for low and intermediate St. We use the model of Zaichik & Alipchenkov (2009) 
to explain the physical mechanisms responsible for these trends, and find that this model 
is able to capture the quantitative behavior of the RDFs extremely well when DNS data 
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for the structure functions are specified, in agreement with Bragg & Collins (2014o). 
We also observe that at large St, changes in the RDF are related to changes the scaling 
exponents of the relative velocity variances. The particle collision kernel closely matches 
that computed by Rosa et al. (2013) and is found to be largely insensitive to the fiow 
Reynolds number. This suggests that relatively low-Reynolds-number simulations may 
be able to capture much of the relevant physics of droplet collisions and growth in the 
adiabatic cores of atmospheric clouds. 


1. Introduction 

Since the pioneering study of Orszag & Patterson (1972o) over forty years ago, direct 
numerical simulation (DNS) has been widely used to study turbulent flows. Previous DNS 
studies have provided a wealth of information about the underlying turbulent fiow field, 
much of which is very difficult to obtain experimentally, including Lagrangian statistics 
(Yeung & Pope 1989), pressure fluctuations (Spalart 1988), and velocity gradient tensors 
(Ashurst et al. 1987). 

Only within the last ten years, however, with the advent of tera- and petascale com¬ 
puting, have DNS at Reynolds numbers comparable to those in the largest laboratory 
experiments become possible. The highest-Reynolds-number simulations to date (with 
Taylor microscale Reynolds numbers R\ ~ 1000) have been of isotropic turbulence in 
tri-periodic domains and have considered both the Eulerian dynamics of the turbulent 
flow field and the Lagrangian dynamics of inertialess tracer (i.e., fluid) particles advected 
by the fiow (Kaneda et al. 2003; Ishihara et al. 2007, 2009; Yeung et al. 2012). 

Many industrial and environmental turbulent flows, however, are laden with dense, 
inertial particles, which can display profoundly different dynamics than inertialess fluid 
particles. The degree to which the dynamics of inertial particles differ from those of fluid 
particles depends on their Stokes number St, a non-dimensional measure of particle iner¬ 
tia, which we define based on Kolmogorov-scale turbulence. We summarize the relevant 
physical mechanisms at small, intermediate, and large values of St below. 

It is well-known from both computational and experimental studies that inertial par¬ 
ticles preferentially sample certain regions of the fiow (e.g., see Balachandar & Eaton 
2010). This preferential sampling is often attributed to the fact that heavy particles are 
centrifuged out of vortex cores and accumulate in low-vorticity and high-strain regions 
(Maxey 1987; Squires & Eaton 1991; Eaton & Fessler 1994), leading to higher colli¬ 
sion rates (Sundaram & Collins 1997). However, this centrifuge mechanism is mainly 
important for small-S't particles which are strongly coupled to the underlying fiow. As 
St is increased, the particle dynamics become less coupled to the local fluid velocity 
field and the influence of their path-history interactions with the turbulence becomes 
increasingly important (e.g., see Bragg & Collins 20146). Particles with sufficiently large 
St can therefore come together from different regions of the fiow with large relative ve¬ 
locities, increasing their collision rate (Wilkinson et al. 2006; Falkovich & Pumir 2007). 
Such a process is referred to as ‘caustics’ (Wilkinson et al. 2006) and the ‘sling effect’ 
(Falkovich & Pumir 2007). At high values of St, several studies (e.g.. Bee et al. 2006o; 
Ayyalasomayajula et al. 2008) have shown that particles have a modulated response to 
the underlying turbulence as they filter out high-frequency fiow features (i.e., features 
with timescales significantly below the particle response time), and they therefore have 
lower kinetic energies and lower accelerations. 

Despite recent advances in simulating high-Reynolds-number turbulent flows, current 
studies of inertial particles in turbulence are primarily at low and moderate Reynolds 
numbers {R\ < 500), and only recently have DNSs been conducted of inertial particles 
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in turbulence with a well-defined inertial range (Bee et al. 2010a,b; Pan et al. 2011; 
Ray & Collins 2011; Rosa et al. 2013; Pan & Padoan 2013). It is vital to understand 
the effect of Reynolds number on the mechanisms above (preferential sampling, path- 
history interactions, and inertial filtering), particularly at higher Reynolds numbers which 
are more representative of those in nature. We give two examples to emphasize the 
importance of developing such an understanding. 

The first example, cloud formation, is the primary motivation for this work. For reviews 
on this subject, see Shaw (2003); Devenish et al. (2012); Grabowski & Wang (2013); 
here we provide a brief overview. It is well-known that standard microphysical cloud 
models over-predict the time required for the onset of precipitation in warm cumulus 
clouds (e.g., see Shaw 2003). At early stages of cloud formation, particles experience 
condensational growth. This process slows down quickly with increasing droplet diameter, 
making condensational growth effective only for droplets with diameters less than about 
OOyim (Grabowski & Wang 2013). Moreover, gravity is only able to significantly enhance 
collisional growth for particles with diameters above 80/im (Pruppacher & Klett 1997; 
Grabowski & Wang 2013), leaving a ‘size gap’ where neither condensational growth nor 
gravitational coalescence is very effective. For particles between these two limits, it has 
been proposed that turbulence-induced collisions are primarily responsible for droplet 
growth. 

It is unclear, however, the extent to which particle collision rates are affected by changes 
in Reynolds number at conditions representative of those in cumulus clouds (which have 
R\ ~ 10,000, see Siebert et al. 2006). Sundaram & Collins (1997) showed that particle 
collision rates depend on both the degree of clustering and on the relative velocities 
between particles, and thus many subsequent analyses have considered the Reynolds- 
number dependence of both of these statistics. While the early study of Wang et al. (2000) 
suggested that clustering increases with R\, later investigations (Collins & Keswani 2004; 
Bee et al. 2010a; Ray & Collins 2011; Rosa et al. 2013) indicate that clustering saturates 
at higher Reynolds numbers. Other researchers have suggested that caustics become more 
prevalent at high Reynolds numbers, leading to larger relative velocities and thus more 
frequent particle collisions (Falkovich et al. 2002; Wilkinson et al. 2006). The findings of 
Bee et al. (2010a) and Rosa et al. (2013), however, do not seem to support that trend. 
In all cases, the Reynolds-number range {R\ < 500) leaves open the question of whether 
the results apply to atmospheric conditions at much higher Reynolds numbers. 

The second example relates to planetesimal formation. Planetesimals begin to form 
when small dust grains collide and coalesce in turbulent protoplanetary nebulae (Pan 
& Padoan 2010). Cuzzi et al. (2001) estimated that the turbulence in such nebulae is 
characterized by R\ ~ 10^ — 10®. It is unclear to what extent the rate of coalescence 
depends on the Reynolds number, and studies at progressively higher Reynolds numbers 
are necessary to develop scaling relations for particle collision rates at conditions repre¬ 
sentative of nebula turbulence. Pan & Padoan (2010) noted that the range of relevant 
particle sizes in the planetesimal formation process spans about nine orders of magnitude, 
and therefore we expect that the collision rates will be affected by preferential sampling 
(for small, medium, and large particles), path-history interactions (for medium and large 
particles), and inertial filtering (for the largest particles). 

In this study, we use high-performance computing resources provided by the U. S. 
National Center for Atmospheric Research (Computational and Information Systems 
Laboratory 2012) to simulate inertial particles in isotropic turbulence over the range 
88 < i?A ^ 597. To our knowledge, the top value represents the highest Reynolds-number 
flow with particles simulated to date. The overall goal is to improve predictions for the 
collision kernel at Reynolds numbers more representative of those in atmospheric clouds. 
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Gravitational forces are neglected in this study, but will be considered in detail in Part 
II of this study (Ireland et al. 2015). 

The paper is organized as follows: §2 provides a summary of the numerical methods 
used and the relevant fluid and particle parameters. In §3, we study single-particle statis¬ 
tics (small-scale velocity gradients, large-scale velocity fluctuations, and accelerations). 
Many of the results from this section help explain the particle-pair statistics presented in 
§4. These statistics include the particle relative velocities, radial distribution functions, 
and collision kernels. Finally, in §5, we summarize our results and suggest practical im¬ 
plications for the turbulence and cloud physics communities. 


2. Overview of simulations 

A brief summary of the simulation parameters and numerical methods is provided 
below. Refer to Ireland et al. (2013) for a more detailed description of the code, including 
integration techniques, parallelization strategies, and interpolation methods. 

2.1. Fluid phase 

We perform DNS of isotropic turbulence on a cubic, tri-periodic domain of length C = 2tt 
with grid points. A pseudospectral method (Orszag & Patterson 1972&) is used to 
evaluate the continuity and momentum equations for an incompressible flow, 

V-M = 0, (2.1) 

^+ujxu + V + =nV^u + f. ( 2 . 2 ) 

dt \pf 2 J 

Here, u is the fluid velocity, a; = V x a is the vorticity, p is the pressure, Pf is the fluid 
density, ly is the kinematic viscosity, and / is a large-scale forcing term that is added 
to make the flow field statistically stationary. For our simulations, we added forcing to 
wavenumbers with magnitude k = in Fourier space in a deterministic fashion to 
compensate precisely for the energy lost to viscous dissipation (Witkowska et al. 1997). 

We perform a series of five different simulations, with Taylor microscale Reynolds 
numbers R\ = 2fc-\/5/ (3:^e) ranging from 88 to 597, where k denotes the turbulent 
kinetic energy and e the turbulent energy dissipation rate. Details of the simulations 
are given in table 1. The simulations are parameterized to have similar large scales, but 
different dissipation (small) scales. The small-scale resolution for the simulations was 
held constant, with Kmax''? ~ 1-6 — 1.7, where Kmax = V^N/3 is the maximum resolved 
wavenumber and rj = {v^/e^ is the Kolmogorov lengthscale. Time-averaged energy and 
dissipation spectra for all five simulations are shown in figure 1. A clear —5/3 spectral 
slope is evident for the three highest Reynolds-number cases {R\ ^ 224), indicating the 
presence of a well-defined inertial subrange. The simulations are performed in parallel on 
Aproc processors, and the P3DFFT library (Pekurovsky 2012) is used for efficient parallel 
computation of three-dimensional fast Fourier transforms. 

2.2. Particle phase 

We simulate the motion of small {d/rj <C 1, where d is the particle diameter), heavy 
{pp/pf ^ 1, where Pp is the particle density), spherical particles. 18 different particle 
classes are simulated with Stokes numbers St ranging from 0 to 30. St = Tp/r,, is a non- 
dimensional measure of a particle’s inertia, comparing the response time of the particle 
Tp = PpdS/ {13pfv) to the Kolmogorov timescale = {vje)^^'^. 
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Table 1. Flow parameters for the DNS study. All dimensional parameters are in arbitrary units, 
and all statistics are averaged over time T. All quantities are defined in the text in §2.1 and 


Simulation 

I 

II 

III 

IV 

V 

R\ 

88 

140 

224 

398 

597 

V 

0.005 

0.002 

0.0008289 

0.0003 

0.00013 

t 

0.270 

0.267 

0.253 

0.223 

0.228 

t 

1.46 

1.41 

1.40 

1.45 

1.43 

(■h 

55.8 

107 

204 

436 

812 

u' 

0.914 

0.914 

0.915 

0.915 

0.915 

u' jUr, 

4.77 

6.01 

7.60 

10.1 

12.4 

Tl 

1.60 

1.54 

1.53 

1.58 

1.57 

Tl/tp 

11.7 

17.7 

26.8 

43.0 

65.4 

T/Tl 

15.0 

10.4 

11.4 

11.1 

5.75 


1.59 

1.59 

1.66 

1.60 

1.70 

N 

128 

256 

512 

1024 

2048 

Np 

262,144 

262,144 

2,097,152 

16,777,216 

134,217,728 

Atracked 

32,768 

32,768 

262,144 

2,097,152 

16,777,216 

proc 

16 

16 

64 

1024 

16,384 




Figure 1. (a) Energy (a) and (b) dissipation spectra for the different simulations described in 
table 1. The diagonal dotted line in (a) has a slope of —5/3, the expected spectral scaling in the 
inertial subrange. All values are in arbitrary units. 


We assume that the particles are subjected to only linear drag forces, which is a reason¬ 
able approximation when the particle Reynolds number Rcp = \u{xP{t),t) — vP{t)\/i' < 
0.5 (Elghobashi & Truesdell 1992). Here, u{x^{f),f) denotes the undisturbed fluid veloc¬ 
ity at the particle position xP{t), and vP{t) denotes the velocity of the particle. (Through¬ 
out this study, we use the superscript p on x, u, and v to denote time-dependent, La- 
grangian variables defined along particle trajectories. Phase-space positions and velocities 
are denoted without the superscript p.) Though particles with large St experience non- 
negligible nonlinear drag forces (e.g., Wang & Maxey 1993), the use of a linear drag 
model for laige-St particles provides a useful first approximation and facilitates compar¬ 
ison between several theoretical models that make the same assumption (e.g., Chun et al. 
2005; Zaichik & Alipchenkov 2009; Gustavsson & Mehlig 2011). The present study also 
neglects the influence of gravity. Part 11 of this study (Ireland et al. 2015) will address the 
combined effects of gravity and turbulence on particle motion. Finally, since a primary 
motivation is to understand droplet dynamics in atmospheric clouds, where the particle 
mass and volume loadings are low (Shaw 2003), we assume that the particle loadings are 
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sufficiently dilute such that inter-particle interactions and two-way coupling between the 
phases are negligible (Elghobashi & Truesdell 1993; Sundaram & Collins 1999). 

Under these assumptions, each inertial particle obeys a simplified Maxey-Riley equa¬ 
tion (Maxey & Riley 1983), 

d'^xP dvP u{xP{t),t) — vP{t) 

^ 

and each fluid (i.e., inertialess) particle is tracked by solving 

f}'r*P 

— =u{xP{t),t). (2.4) 

To compute uP{t) = u{xP{t),t), we need to interpolate from the Eulerian grid to the 
particle location. While other studies (e.g., see Bee et al. 2010a; Durham et al. 2013) have 
done so using tri-linear interpolation, Ireland et al. (2013) showed that such an approach 
can lead to errors in the interpolated velocity which are orders of magnitude above the 
local time-stepping error. In addition, van Hinsberg et al. (2013) demonstrated that tri- 
linear interpolation, which possesses only C° continuity, leads to artificial high frequency 
oscillations in the computed particle accelerations. Ray & Collins (2013) noted that the 
relative motion of particles at small separations will depend strongly on the interpolation 
scheme. Since a main focus of this paper is particle motion near-contact and its influence 
on particle collisions, it is crucial to calculate u{xP{t),t) as accurately as possible. To 
that end, we use an eight-point B-spline interpolation scheme (with C® continuity) based 
on the algorithm in van Hinsberg et al. (2012). 

The particles were initially placed in the flow with a uniform distribution and velocities 
vP equal to the underlying fluid velocity uP. We began computing particle statistics once 
the particle distributions and velocities became statistically stationary, usually about 5 
large-eddy turnover times = £/u' (where £ is the integral lengthscale and u' = y/2k/3) 
after the particles were introduced into the flow. Particle statistics were calculated at a 
frequency of 2-3 times per and were time-averaged over the duration of the run T. 

For a subset A^tracked of the total number of particles in each class Np, we stored particle 
positions, velocities, and velocity gradients every O.It^ for a duration of about lOOr,,. 
These data are used to compute Lagrangian correlations, accelerations, and timescales 
of the particles. 


3. Single-particle statistics 

We first consider single-particle statistics from our simulations. These statistics will 
provide a basis for our understanding of the two-particle statistics presented in §4. We 
explore velocity gradient (i.e., small-scale velocity) statistics in §3.1, kinetic energy (i.e., 
large-scale velocity) statistics in §3.2, and acceleration statistics in §3.3. In each case, we 
study the effect of the underlying flow topology on these statistics. 

3.1. Velocity gradient statistics 

We consider the gradients of the underlying fluid velocity at the particle locations, 
A(xP{t),t) = Wu{xP{t),t). These statistics provide us with information about the small- 
scale velocity field experienced by the particles. (Refer to Meneveau (2011) for a recent 
review on this subject.) In particular, to understand the interaction of particles with 
specific topological features of the turbulence, we decompose A{xP{t),t) into a symmet¬ 
ric strain rate tensor S{xP{t),t) = [A{xP{t),t) + A'^ {xP (t), t)]/ 2 and an antisymmetric 
rotation rate tensor 7t{xP{t),t) = [A{xP{t),t) — A'^(xP(t),t)]/‘2. 
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Due to their inertia, heavy particles are ejected out of regions of high rotation rate 
and accumulate in regions of high strain rate (e.g., Maxey 1987; Squires & Eaton 1991; 
Eaton & Fessler 1994), and this is associated with a ‘preferential sampling’ of A{x,t). 
For particles with low inertia {St 1), preferential sampling is the dominant mechanism 
affecting the particle motion (e.g., see Chun et al. 2005). As the particle inertia increases, 
the particle motion becomes increasingly decoupled from the local fluid turbulence, and 
the effect of the preferential sampling on the particle dynamics decreases. At the other 
limit {St ^ 1), preferential sampling vanishes and the particles have a damped response 
to the underlying flow which leads them to sample the turbulence more uniformly (e.g., 
see Bee et al. 2006a). 

We first consider the average of the second invariants of the strain rate and rotation 
rate tensors evaluated at the inertial particle positions 


{Sy = {S{xP{t),t): S{xP{t),t)), (3.1) 

and 

(7^2)P = {n{xP{t),t): n{xP{t), t)). (3.2) 

By definition, for fully mixed fluid particles {St = 0) in homogeneous turbulence, t{^{S‘^)p = 
T^{n^)p = 0.5. 

Since small-S't particles are centrifuged out of regions of high rotation, we expect that 
T{^{Vf)P will decrease with increasing St; their accumulation in high strain regions would 
also lead to the expectation that t‘^{S‘^)p will increase with increasing St. In figure 2 we 
see that while t{^{TZ‘^)p is more strongly affected by changes in Rx than is r^(5^)^, both 
quantities decrease with increasing St (for St 1). This surprising result is consistent 
with other DNS (Collins & Keswani 2004; Chun et al. 2005; Salazar & Collins 2012o). 
Our data also show that both r^(5^)^ and t{^{TZ‘^)p decrease with increasing Rx for 
St 1, in agreement with Collins & Keswani (2004). 

We use the formulation given in Chun et al. (2005) (and re-derived in Salazar & 
Collins 2012a) to model the effect of preferential sampling on t‘^{S‘^Y and t{^{TZ'^)p in 
limit of St 1. Chun et al. (2005); Salazar & Collins (2012a) showed that for an 
arbitrary quantity the average value of (j) sampled along a particle trajectory {(j))P can 
be reconstructed entirely from fluid particle statistics using the relation, 

{(j){St)Y = {(j){St = 0))P -f T^crlSt • (3.3) 

Here, Uy denotes the standard deviation of a variable Y along a fluid particle trajectory, 
Pyz is the correlation coefficient between Y and Z, 


P _ 

Pyz = 


[Y{xP{t),t) - {Y{xP{t),t))] [Z{xP{f),t) - {Z{xP{t),t))] 


G\rG 


(3.4) 


Y'-' Z 


and Ty^ is the Lagrangian correlation time. 


rpP _ J 0 
^YZ — 


[Y{x^{0),0) - {Y{x^{t),t))] [Z{xP{t'),t') - {Z{xP{t),t))]) dt' 


\Y{xP{t),t) - {Y{xP{t),t))] [Z{xP{t),t) - {Z{xP{t),t)) 


(3.5) 


The predictions from (3.3) for small St are shown by the solid lines in figure 2(c) and 
figure 2(d). In the limit of small St, this model is able to capture the decrease in both 
t{^{S‘^Y and t{^{TYY with increasing St, and also the decrease in these quantities with 
increasing Rx. It is uncertain whether the quantitative differences between the DNS data 



P. J. Ireland et al. 



St 




St 



St 


Figure 2. Data for (a,c) and {TZ^)^ (b,d) sampled at inertial particle positions as function 

of St for different values of R\. The data are shown at low St in (c,d) to highlight the effect of 
preferential sampling in this regime. The solid lines in (c) and (d) are the predictions from (3.3) 
for St 1. DNS data are shown with symbols. 


and the model are due to shortcomings of the model or the fact that the smallest inertial 
particles {St = 0.05) are too large for the model (which assumes St <C 1) to hold. 

Despite the success of the model of Chun et al. (2005) in reproducing the trends 
in the DNS, the physical explanation for the changes in the mean strain and rotation 
rates remains unclear. In figure 3(a), we plot joint PDFs of the strain and rotation 
rates sampled by both St = 0 and St = 0.1 particles to better understand the specific 
topological features of the regions of the flow contributing to these changes. Following 
the designations given in Soria et al. (1994), we refer to regions with high strain and 
high rotation (indicated by ‘A^ in figure 3(a)) as ‘vortex sheets,’ regions of low rotation 
and high strain (indicated by ‘B’) as ‘irrotational dissipation’ areas, and regions of high 
rotation and low strain (indicated by ‘C”) as ‘vortex tubes.’ 

Our results show three main trends in the particle concentrations. First, inertial par¬ 
ticles are ejected from vortex sheets (A) into regions of moderate rotation and moderate 
strain {A'). This ejection from vortex sheets has only recently been discussed in the liter¬ 
ature (Salazar & Collins 2012a). Second, they move from irrotational dissipation regions 
[B] to regions of comparable rotation and even higher strain {B'). Third, the particles 
move out of vortex tubes (C) into regions of lower rotation and higher strain {C). Ev¬ 
idently, this first effect is primarily responsible for the decrease in small St, 

as suggested in Salazar & Collins (2012a), and the first and third effects both contribute 
to the decrease in We will revisit these three trends in relation to the particle 

kinetic energies (§3.2) and the particle accelerations (§3.3). 
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Figure 3. (a) Joint PDFs of r^<S: S and 'R, for R\ = 597 for St = 0 and St = 0.1 

particles. Certain regions of the flow are labeled to aid in the discussion of the trends, (b) Joint 
PDFs of : S and TZ. for different Rx for St = 0 particles. In both plots, the exponents 

of the decade are indicated on the contour lines. 



Figure 4. The difference between the mean rates of strain and rotation sampled by the 
particles as a function of St for different values of Rx- 

Figure 3(b) shows the PDF map for fluid particles at three values of the Reynolds 
number. Notice that as Rx increases, the probability of encountering a vortex sheet 
(overlapping high strain and high rotation) increases. This finding is consistent with the 
results of Yeung et al. (2012), who observed that high strain and rotation events increas¬ 
ingly overlap in isotropic turbulence as the Reynolds number increases. It is thus likely 
that with increasing Reynolds number, rotation and strain events become increasingly in¬ 
tense, and the resulting vortex sheets become increasingly efficient at expelling particles, 
causing both and to decrease (cf. figure 2). 

Maxey (1987) noted that at low St, the compressibility of the particle field (and hence 
the degree of particle clustering) is directly related to the difference between the rates 
of strain and rotation sampled by the particles, T^{S^y — r^(7^^)^. From figure 4, we 
see that at low St, increases with increasing Rx, suggesting that the 

degree of clustering may also increase here. We will test this hypothesis in §4.2 when we 
directly measure particle clustering at different values of St and Rx- 
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Figure 5. Lagrangian timescales of a single component of the strain rate (a) and rotation rate 
(b) tensors, plotted as a function of St for different values of Rx- 


We finally consider the Lagrangian strain and rotation timescales, which will be useful 
for understanding the trends in particle clustering in §4.2. Since the fluid and particle 
phases are isotropic, we will have nine statistically equivalent strain timescales: ^ , 

s jTS o , To o , To o , To o , To o , To o , and To o . We take the strain 
timescale to be the average of these nine components. We similarly take the rotation 
timescale to be the average of three statistically equivalent components: T^ ^ , 

We see that is independent of R\ for St < 10, and decreases weakly with 

increasing Rx for St > 10. On the other hand, T^^/rn tends to decrease with increasing 
Rx for all values of St, and this decrease becomes more pronounced as St increases. We 
also see that T^.^^ is much more sensitive to changes in St than Tgg, suggesting that 
the dominant effect of inertia is to cause particles to spend less time in strongly rotating 
regions. As a result, the particles will generally have less time to respond to fluctuations 
in the rotation rate, causing to be strongly reduced with increasing St, as was seen 

above. 


3.2. Particle kinetic energy 

We now move from small-scale velocity statistics to large-scale velocity statistics. Figure 6 
shows the average particle kinetic energy kP{St) = ^{vP{t) ■ v^{t)) (normalized by the 
average fluid kinetic energy k) for different values of Rx- 

We first consider the effect of inertial filtering on this statistic, and then examine the 
effect of preferential sampling. It is well-known that filtering leads to a reduction in the 
particle turbulent kinetic energy for large values of St. This reduction is the strongest 
(weakest) for the lowest (highest) Reynolds numbers, as seen in figure 6(a). These trends 
are captured by the model in Abrahamson (1975), which assumes an exponential decor¬ 
relation of the Lagrangian fluid velocity. Under this assumption, the ratio between the 
particle and fluid kinetic energies can be expressed as 

k^St) ^ 1 ^ 1 ,3 

k l + Tp/n 1 + St{Tri/n)’ 

where is the Lagrangian correlation time of the fluid, which we approximate using the 
relation given in Zaichik et al. (2003). The model predictions of k^{St)/k are included 
in figure 6(a) and are in good agreement with the DNS at large St, where filtering is 
dominant. The trends with Rx are also reproduced well. 

We thus have the following physical explanation of inertial filtering on the particle 
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Figure 6. (a) The ratio between the average particle kinetic energy k^{St) and the average 
fluid kinetic energy k for different values of R\. DNS data are shown with symbols, and the 
predictions of the filtering model in (3.6) are shown with solid lines, (b) The ratio between 
and k (open symbols), and the ratio between the average fluid kinetic energy at the 
particle locations k^^{St) and k (filled symbols), shown at low St to highlight the effects of 
preferential sampling. Also shown is the prediction from the preferential sampling model given 
in (3.3) (solid lines). 

kinetic energies: for low-Reynolds-number flows, the response time of the largest particles 
exceeds the timescales of many large-scale flow features. The result is a Altered response 
to the large-scale turbulence and an overall reduction in the particle kinetic energy. As 
the Reynolds number is increased (and the particle response time is fixed with respect to 
the small-scale turbulence), more flow features are present with timescales that exceed 
the particle response time, and hence the effect of inertial filtering is diminished with 
increasing R\, as predicted by (3.6). 

To highlight the effect of preferential sampling on the particle kinetic energy, fig¬ 
ure 6 (b) shows both the average particle kinetic energy k^^St) and the average kinetic 
energy of the fluid sampled along an inertial particle trajectory, yP{Si) = ^{u{xP{t),t) ■ 
u{xP{t),t)). As is evident in figure 6 (b), the particle kinetic energy exceeds k for low 
values of St. By comparing k^ to k^P, we see that the increased kinetic energy of the 
smallest particles is due almost entirely to preferential sampling of the flow field. While 
Salazar & Collins (20126) were the first to show an increase in kP{St)/k for low St (which 
they attributed to preferential sampling), this trend is also suggested by the early study 
of Squires & Eaton (1991), in which the authors observed that small inertial particles 
preferentially sample certain high kinetic energy regions they referred to as ‘streaming 
zones.’ Figure 6 (b) also shows that at small values of St, kP{Sf)/k decreases with in¬ 
creasing Reynolds number. 

The solid lines in figure 6 (b) show the predictions of the particle kinetic energy from 
(3.3). In the limit of small St, the model of Chun et al. (2005) is able to capture quali¬ 
tatively both the increase in kP{St)/k with increasing St and the decrease in kP{St)/k 
with increasing R\. 

To further elucidate the physical mechanisms leading to these trends, we plot the mean 
kinetic energy of the fluid conditioned on and , ^52 7 ^ 2 , in figure 7. Isocontours of 
the concentrations oi St = 0 and St = 0.1 particles are shown for comparison. While 
the data contain considerable statistical noise, we can draw a few conclusions about the 
qualitative trends. 

From figure 7(a), we see that the change in kinetic energy at R\ = 88 can be divided 
into the three mechanisms discussed in §3.1. First, particles are ejected from vortex sheets 
(A) into moderate rotation and moderate strain regions {A'), which generally tends to 
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Figure 7. Filled contours of the fluid kinetic energy conditioned on and TZ^, ^52 7 ^ 2 , normal¬ 
ized by the unconditioned mean fluid kinetic energy k at (a) the lowest Reynolds number and 
(b) the highest Reynolds number. The dotted contour lines indicate kg 2 -jz^/k = 1. Isocontours 
of particle concentration for St — 0 and St = 0.1 particles are included for reference, with the 
exponents of the decade indicated on the contour lines. Certain regions of the flow are labeled 
to aid in the discussion of the trends. 

decrease the particle kinetic energy. Second, as St increases, particles in irrotational 
straining regions (B) travel into regions of higher strain {B'), which are characterized 
by higher kinetic energy. Third, some inertial particles are ejected from vortex tubes 
(C), which are characterized by lower kinetic energies, and travel into lower rotation and 
higher strain regions (C"), which have higher kinetic energies. The observed increase in 
kP{St)/k must therefore be due to the second and third mechanisms. 

At high Reynolds numbers (figure 7(b)), however, a larger portion of the flow is oc¬ 
cupied by regions of overlapping high strain and high rotation from which particles are 
ejected (see §3.1). The first mechanism (which tends to decrease the kinetic energy) 
therefore plays a larger role. Also, at R\ = 597, high rotation and low strain regions (C) 
are no longer associated with very low kinetic energies, causing the third mechanism to 
be less effective at increasing the particle kinetic energy. The overall result is a decrease 
in kP{St)/k with increasing Reynolds number at small values of St. 

3.3. Particle accelerations 

In this section, we analyze fluid and inertial particle accelerations a^(t) = dvP{t)/dt. 
Fluid particle accelerations are known to be strongly intermittent (e.g., see Voth et al. 
2002; Ishihara et al. 2007), with the probability of intense acceleration events increasing 
with the Reynolds number. Before accounting for inertial effects, we consider the effect 
of R\ on the acceleration variance (a^)^ = {aP{t) • a‘’{t))/3 of Lagrangian fluid particles 
in figure 8(a). To facilitate comparison between the different Reynolds numbers, we have 
normalized (a^)^ by the Kolmogorov acceleration variance The DNS data 

from Yeung et al. (2006) and the theoretical predictions of Hill (2002), Sawford et al. 
(2003), and Zaichik et al. (2003) are shown for comparison. We see that our DNS data 
agrees well with Yeung et al. (2006), and that the model of Sawford et al. (2003) best 
reproduces the trends in the DNS. Hill (2002) breaks down at low R\, while Zaichik et al. 
(2003) fails at high R\. 

We turn our attention to inertial particle accelerations in figure 8(b). The observed 
trend for inertial particles is analogous to that for fluid particles: at each value of St 
considered, the particle acceleration variance (normalized by Kolmogorov units) mono- 
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Figure 8. (a) The acceleration variance of Lagrangian fluid particles as a function of R\. The 
results from the present study (open circles) are compared to DNS data from Yeung et al. 
(2006) (filled squares) and several theoretical predictions (lines), (b) The acceleration variance 
of inertial particles as a function of St for different values of R\. 


tonically increases with R\ (cf. Bee et al. 2006a). As St increases, the acceleration vari¬ 
ance decreases, presumably as a result of both preferential sampling of the flow field and 
inertial filtering. 

We now seek to understand and model how inertia changes the accelerations of particles 
through the hltering and preferential sampling effects. To do so, we rescale the inertial 
particle acceleration variance by that of fluid particles and plot the results in figure 9. 
In figure 9(a), we compare the rescaled acceleration variance to the model of Zaichik & 
Alipchenkov (2008), which only accounts for inertial filtering of the underlying flow. The 
model of Zaichik & Alipchenkov (2008) is able to capture all the qualitative trends in R\ 
and St, and the model predictions provide remarkably good quantitative agreement with 
the DNS at the largest values of St, where filtering is the dominant mechanism. At lower 
values of St, the rescaled particle acceleration variance decreases with increasing Rx- In 
this case, as Rx increases, the underlying flow is subjected to increasingly intermittent 
acceleration events, and the inertial particles filter a larger fraction of these events. At the 
largest values of St, most intermittent accelerations are hltered, and a particle’s acceler¬ 
ation variance is determined by its interaction with the largest turbulence scales. Since 
the range of available large scales increases with Rx, the rescaled particle acceleration 
variance increases with Rx for the largest values of St. 

We now consider the effect of preferential sampling on the acceleration variances. In 
figure 9(b), we plot the variance of both inertial particle accelerations and fluid acceler¬ 
ations along inertial particle trajectories (scaled by the acceleration variance of S't = 0 
particles). As expected, for S't <C I, where preferential sampling is the dominant mech¬ 
anism, inertial particle accelerations are almost equivalent to the accelerations of the 
underlying flow sampled along the particle trajectories. The model of Chun et al. (2005) 
(3.3) is able to reproduce all the qualitative trends correctly in the limit of small St. The 
scaled variances decrease with increasing Rx, and we expect that this trend is due to 
the fact that high vorticity regions are associated with high accelerations (Biferale et al. 
2005) and become increasingly efficient at ejecting particles (refer to §3.1). 

We test this expectation in figure 10 by plotting the acceleration variance for fluid 
particles conditioned on and Tif, (a ^)^2 7 ^ 2 , and normalized by the unconditioned 
variance {af)^. We see that inertial particles are indeed ejected from high vorticity re¬ 
gions (both vortex sheets and vortex tubes) into lower vorticity regions (e.g., A into A' 
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Figure 9. (a) Inertial particle acceleration variances scaled by the fluid particle acceleration 
variance (open symbols). The solid lines and arrows indicate the predictions from the filtering 
model of Zaichik & Alipchenkov (2008). (b) The variance of the inertial particle accelerations 
(open symbols) and the flnid velocity accelerations along the particle trajectories (filled symbols), 
shown at low St to highlight the effect of preferential sampling. The solid lines indicate the 
predictions from the preferential sampling model given in (3.3). 



K:K K:K 

Figure 10. Filled contours of the variance of the fluid particle accelerations conditioned on 
and TZ^, (a ^)^2 ^^ 2 , normalized by the unconditioned fluid particle acceleration variance {a?)^, 
at (a) R\ = 88 and (b) R\ = 597. Isocontours of particle concentration for St = 0 and St = 0.1 
particles are included for reference, with the exponents of the decade indicated on the contour 
lines. Certain regions of the flow are labeled to aid in the discussion of the trends. 


and C into C"), and that these high vorticity regions are marked by very large accel¬ 
erations. Though some inertial particles experience higher accelerations as they move 
into irrotational straining regions with higher strain rates (e.g., B into B'), this effect 
is relatively weak, and the overall trend is a decrease in the particle accelerations with 
increasing inertia. 

To investigate the intermittency of inertial particle accelerations, we plot the kurtosis 
of the particle accelerations, {a'^Y/{{a?YY, in figure 11, where (a^)^ = {a\{tY + a^YY + 
o^YYY/3. (Note that a Gaussian distribution has a kurtosis of 3, as indicated in figure 11 
by a dotted line.) As expected, the particle accelerations are highly intermittent, with 
the degree of intermittency increasing with increasing R\. The kurtosis decreases very 
rapidly as St increases. Figure 11(b) indicates that the kurtosis of very small particles 
{St = 0.05) at the highest value of R\ is over a factor of two smaller than that of fluid 
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Figure 11. Particle acceleration kurtosis as a function of St for different values of Rx. The 
dotted line indicates a kurtosis of 3, the value for a Gaussian distribution. Values over the 
whole range of non-zero St are shown in (a), (b) shows only small-St results on a linear plot to 
emphasize the rapid reduction in kurtosis as St increases from 0. 


particles. The largest-S”! particles have kurtosis values approaching those of a Gaussian 
distribution. These trends can be explained by the fact that both preferential sampling 
and inertial hltering decrease the probability of high-intensity acceleration events. Stan¬ 
dardized moments of up to order 10 (not shown) were also analyzed and found to exhibit 
the same trends. 

We should note that the grid resolution study in Yeung et al. (2006) suggests that 
the acceleration moments from our DNS may be under-predicted. Yeung et al. (2006) 
showed that at Rx « 140, increasing the grid resolution fcmax?? from 1.5 to 12 led to a 10% 
increase in the fluid acceleration variance and a 30% increase in the fluid acceleration 
kurtosis. It is unclear how these trends will change at higher Rx, but it suggests that 
the quantitative results reported here should be interpreted with caution. (The velocity 
gradients presented earlier are likely reliable, however, since Yeung et al. (2006) found 
that such statistics are less dependent on the grid resolution.) 


4. Two-particle statistics 

We now consider two-particle statistics relevant for predicting inertial particle colli¬ 
sions. We analyze particle relative velocities in §4.1, clustering in §4.2, and use these data 
to compute the collision kernel in §4.3. (The mean-squared separation of inertial particle 
pairs was also studied from these data and is the topic of a separate publication (Bragg 
et al. 2015a).) 

4.1. Particle relative velocities 

We study particle relative velocities as a function of both St and Rx. The relative veloc¬ 
ities for inertial particles are defined by the relation 

wl^{t)^[vl{t)-vl{t)]-el^{t). (4.1) 

Here, v\ and uf indicate the velocities of particles 1 and 2, respectively, which are sep¬ 
arated from each other by a distance r^{f) = \rt’{t)\. The subscripts jj and T indicate 
directions parallel (longitudinal) to the separation vector or perpendicular (transverse) 
to the separation vector, respectively, and denotes the unit vector in the correspond¬ 
ing direction. (We use the method discussed in Pan & Padoan (2013) to compute the 
transverse components.) 
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We will also examine the velocity differences of the fluid at the particle locations, 
defined as 

Au||__L(t) = [uPit) - uP{t)] ■ _L(t), (4.2) 

where and are the velocities of the fluid underlying particles 1 and 2, respectively. 
Note that for uniformly-distributed fluid {St = 0) particles, the particle velocity statistics 
are equivalent to the underlying fluid velocity statistics. 

Following the nomenclature in Bragg & Collins (2014a,&), we denote particle relative 
velocity moments of order n as 

for the components parallel to the separation vector, and as 

(4.4) 

for components perpendicular to the separation vector. In these expressions (•)r denotes 
an ensemble average conditioned on r^(t) = r. 

For the purposes of computing the collision kernel (see §4.3), we are also interested in 
the mean inward relative velocity parallel to the separation vector, defined as 

-5'0||(r) = - / w||p(w|||r)dw||, (4.5) 

J —OO 

where p('u;|||r) = (5('u;||'(t) —W||))r is the PDF for the longitudinal particle relative velocity 
conditioned on r^’(t) = r. 

Finally, in some cases we are also interested in moments of the fluid velocity differences. 
We use a superscript fp to denote the moments of fluid velocity differences at the particle 
locations, and a superscript / to denote the moments of fluid velocity differences at fixed 
points with separation r. We therefore have 

= ( [^W||(rP(t),t)]”^^, (4.6) 

and 

= ( [AM||(r,t)]”). (4.7) 

The components perpendicular to the separation vector are defined analogously. 

We consider dissipation-range statistics in §4.1.1 and inertial-range statistics in §4.1.2. 

4.1.1. Dissipation range relative velocity statistics 

In figure 12, we plot the relative velocity variances and versus r/p at R\ = 597. 
The mean inward relative velocity (not shown) has the same qualitative trends, and will 
be considered later in this section. For the purposes of the following discussion, we define 
the dissipation range as the region over which the fluid velocity variances follow re¬ 
scaling, which is seen to be 0 ^ r/p < 10 in figure 12, in agreement with Ishihara et al. 
(2009). 

At small separations, the relative velocity variances parallel to the separation vector 
(figure 12(a)) increase monotonically with St and deviate from re-scaling, while the 
relative velocity variances perpendicular to the separation vector decrease for At 0.1 
and then increase monotonically with St for St > 0.1 (figure 12(b)). We expect that the 
trends at small separations and small St are primarily due to preferential sampling of 
the underlying flow, which also dictates much of the single-particle dynamics for small 
St (refer to §3). 
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Figure 12. The particle relative velocity variances parallel to the separation vector (a) and 
perpendicular to the separation vector (b), plotted as a function of the separation r/y for 
R\ — 597. The Stokes numbers are indicated by the line labels, and the St = 0 curves are 
shown with dashed lines for clarity. The expected dissipation and inertial range scalings (based 
on Kolmogorov 1941) are included for reference. 




Figure 13. The parallel (a) and perpendicular (b) relative velocity variances of inertial particles 
(Sjii and open symbols) and of the fluid at inertial particle positions {S^^ and filled 
symbols) for R\ = 597. All quantities are normalized by the relative velocity variances of St = 0 
particles. 


To test this expectation, we compare the particle relative velocity variances to those 
of the fluid sampled by the particles in figure 13. In all cases, the velocity variances 
are normalized by those of S't = 0 particles. At St = 0.05 and St = 0.1, the effect of 
preferential sampling is dominant at all separations, as evidenced by the fact that 

and S'2_5^ are close to and 5'2_|_, respectively. We note that for small St and small r/rj, 

preferential sampling leads to an increase in with increasing St and to a decrease 

in Slf_ with increasing St. This is consistent with the trends observed in figure 12 and 
with our argument (§3.1) that inertia causes particles to be ejected from vortex tubes. 
We expect that two particles which are rotating in a vortex tube will experience small 
(large) relative velocities parallel (perpendicular) to the particle separation vector, and 
that the parallel (perpendicular) relative velocities will increase (decrease) as particles 
are ejected from a vortex tube. 

For St > 0.2, the particle relative velocities are much larger than the underlying fluid 
velocity differences at small separations. This difference is due to path-history effects 
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(see Bragg & Collins 2014a,6). That is, as inertial particles approach each other, they 
retain a memory of more energetic turbulence scales along their path histories, leading 
to relative velocities that exceed the local fluid velocity difference. These path-history 
effects imply that inertial particles can come together from different regions in the flow, 
occupy the same position in the flow at the same time, and yet have different velocities 
due to their differing path histories. This effect is referred to as ‘caustics,’ ‘crossing 
trajectories,’ or ‘the sling effect,’ causes a departure from r^-scaling in the second-order 
structure functions at small separations, and can lead to large relative velocities (Yudine 
1959; Falkovich et al. 2002; Wilkinson & Mehlig 2005; Wilkinson et al. 2006; Falkovich 
& Pumir 2007). (Also note that while caustics are instantaneous events, the statistical 
manifestation of caustics is known as ‘random, uncorrelated motion’ and is discussed in 
IJzermans et al. (2010).) Since the timescale over which the particles retain a memory 
of their interactions with turbulence increases with increasing inertia, caustics become 
more prevalent as St increases. 

One effect of caustics is to make the parallel and perpendicular relative velocity com¬ 
ponents nearly the same in the dissipation range, as can be seen in figure 12 for St > 0.3. 
(Note that fluid particles do not experience caustics and have 2S'2|| = S'^.l ^ 1 

as a result of continuity (e.g., see Pope 2000).) For St ^ 10, the relative velocities are 
almost unaffected by the underlying turbulence in the dissipation range. As a result, the 
relative velocities are nearly independent of r/rj in this range. 

The effect of caustics can also be clearly seen in figure 14(a,b), where we plot the 
parallel relative velocities at a given separation as a function of St. From this figure, it 
is evident that the particle relative velocities at the smallest separation sharply increase 
as St exceeds about 0.2. The rapid increase in the particle relative velocities with St is 
consistent with the notion that caustics take an activated form (Wilkinson et al. 2006) and 
that they are negligible below a critical value of St (Salazar & Collins 20126; IJzermans 
et al. 2010). Our data suggest a critical Stokes number for caustics of about 0.2 to 0.3, in 
agreement with Falkovich & Pumir (2007) and Salazar & Collins (20126). The increase 
in the relative velocities occurs at higher values of St as the separation increases. In this 
case, the particles are subjected to larger-scale turbulence, and hence the particles must 
have more inertia for their motion to deviate significantly from that of the underlying 
flow. 

We now examine the Reynolds-number dependence of the relative velocities, restricting 
our attention to the component parallel to the separation vector. The relative velocities 
of the largest particles {St > 10) increase strongly with increasing R\ in figure 14(a,b). 
There are two reasons for this trend. The first is that the effect of filtering on the larger 
turbulence scales decreases as Rx is increased (see §3.2). The second is that u'/u,^ in¬ 
creases with increasing Rx, indicating that large-St particles in the dissipation range 
carry a memory of increasingly energetic turbulence (relative to the Kolmogorov scales) 
in their path history as Rx is increased. 

For smaller values of St {St < 3), the relative velocities in figure 14(a,b) are only 
weakly dependent on Rx, in agreement with previous DNS studies (Wang et al. 2000; 
Bee et al. 2010a; Rosa et al. 2013; Onishi et al. 2013; Onishi & Vassilicos 2014) and the 
model of Pan & Padoan (2010). To highlight any small Reynolds-number effects in this 
range, we therefore divide the relative velocities at r = 0.25?7 and a certain Rx by their 
value at Rx = 88 and plot the results in figure 14(c,d). 

For St < 1, the relative velocity variances increase weakly with increasing Rx (fig¬ 
ure 14(d)). As Rx increases, the range of velocity scales u'/ur^ increases, allowing some 
particle pairs in this Stokes-number range to sample more energetic turbulence as they 
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Figure 14. (a) The mean inward relative velocities and (b) the relative velocity variances, 
plotted as a function of St for small separations and different values of R\. Open symbols denote 
r = 0.25y, gray filled symbols denote r = 1.75?7, and black filled symbols denote r = Q-Tby. To 
emphasize any Reynolds-number dependencies for St ^ 3, we also plot (c,d) the ratio between 
the value of these quantities at a given Reynolds number to their value at Rx — 88 at separation 
r = 0.25?7. 

converge to small separations, and causing the the relative velocity variances to increase 
with increasing Reynolds number. Furthermore, turbulence intermittency, which also 
increases with increasing Reynolds number, may also contribute to the trend in the rel¬ 
ative velocity statistics. We note that the mean inward velocities (figure 14(c)) are less 
affected by changes in Reynolds number, presumably because the mean inward velocity 
is a lower-order statistic that is less influenced by the relatively rare events described 
above. 

For 1 ^ S'! < 3, we also expect the increased scale separation, the increased intermit¬ 
tency of the turbulence, or both to act to increase the relative velocities. However, we 
observe an overall decrease in the relative velocities with increasing R\ here, in agreement 
with Bee et al. (2010a); Rosa et al. (2013). These reduced relative velocities are likely 
linked to the decrease in the Lagrangian rotation timescales with increasing R\ 

observed in §3.1. That is, as decreases with increasing Rx, the particles have 

a shorter memory of fluid velocity differences along their path histories, which in turn 
causes the relative velocities to decrease. _ ^ 

We now examine the behavior of the scaling exponents of oc r^n and oc r^n 
at small separations. (These scaling exponents will also be used in §4.2 to understand 
and predict the trends in the particle clustering.) We compute and Cy using a linear 
least-squares regression for 0.75 < r/rj < 2.75 at different values of St and Rx. Note that 









20 


P. J. Ireland et al. 



Figure 15. Dissipation-range scaling exponents for (a) and for various values of St 
and R\. The exponents are computed from linear least-squares regression for 

0.75 s; r/ri ^ 2.75. 


while using such a large range of r /77 will necessarily introduce finite-separation effects, 
there is generally too much noise in the data to accurately compute the scaling exponents 
over smaller separations. 

The scaling exponents are plotted in figure 15. We note that the scaling exponents are 
below those predicted by Kolmogorov (1941) (hereafter ‘K4T) for fluid {St = 0) particles 
(Cl)" = 1 and Cy = 2) and, like the relative velocities themselves, vary only slightly as R\ 
changes. 

For St ^ 10, the scaling exponents are about zero, indicating that the relative velocities 
are generally independent of r, as explained above. The scaling exponents for 1 ^ S’! 3 

generally increase with increasing R\, since path-history interactions (which generally 
decrease the scaling exponents) become less important, as explained above. Finally, we 
note that Cy decreases with increasing R\ for St ^ 1, since intermittent path-history 
effects are expected to be more important here. 

We next consider the PDFs of the relative velocities in the dissipation range. Figure 16 
shows the PDFs for 0 < r /?7 < 2 and R\ = 597. In figure 16(a), we see that as St increases, 
the tails of the PDF of rcy’/rt,) become more pronounced, indicating that larger relative 
velocities become more frequent, in agreement with our observations above. 

We show PDFs in standardized form in figure 16(b) to analyze the extent to which 
they deviate from that of a Gaussian distribution. It is evident that the degree of non- 
Gaussianity peaks for St ^ 1 and becomes smaller as St increases. The physical expla¬ 
nation for this intermittency at S'! ~ 1 is that the motion of these particles is affected by 
both the small-scale underlying turbulence and by the particles’ memory of large-scale 
turbulent events in their path histories. This combination of contributions from both 
large- and small-scale events leads to strong intermittency. We also see that the underly¬ 
ing fluid is itself quite intermittent at this small separation, as expected (e.g., see Gotoh 
et al. 2002 ). 

We now use three statistical measures to quantify the shape of the PDFs. The first 
is the ratio between the mean inward relative velocities and the standard deviation of 
the relative velocities, *S'((y/( 5 ' 2 y)^/^; the second is the skewness of the relative velocities, 
third is the kurtosis of the relative velocities, 'S' 4 y/(*S' 2 ||)^' (Due 
to insufficient statistics, we will not consider data from these latter two quantities for 
r/rj < 1.75.) 

We show the ratio 5 '((||/(S' 2 y)^/^ in figure 17. One motivation for looking at this ratio 
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Figure 16. PDFs of the particle relative velocities for separations 0 ^ r/y ^ 2 and Rx = 597. 
The relative velocities are normalized by both (a) and (b). The solid lines denote 

the relative velocity PDFs for St — 0 particles, and the dotted line in (b) indicates a standard 
normal distribution. 



Figure 17. The ratio between mean inward relative velocities and the standard deviation of 
the relative velocities as a function of St for small separations and different values of Rx- Open 
symbols denote r = 0.2577, gray filled symbols denote r — 1.7577, and black filled symbols 
denote r = 9.75?;. The horizontal dotted line indicates that value of this quantity for a Gaussian 
distribution. 

is that existing theories (e.g., see Zaichik et al. 2003; Pan & Padoan 2010) only predict 
the relative velocity variance, and by assuming the relative velocities have a Gaussian 
distribution, relate this variance to the mean inward relative velocity. For a Gaussian 
distribution, this ratio is approximately 0.4. At all values of St, R\, and r/rj, our data 
indicate that the ratio is below 0.4 and thus that the particle relative velocities are 
intermittent (see also Wang et al. 2000; Pan & Padoan 2013). The degree of intermittency 
peaks for order unity St, high Rx, and small r/rj, and using a Gaussian prediction in this 
regime would lead to predictions of the mean inward velocity which are in error by more 
than a factor of 2. 

We next consider the skewness, «5'f||/(-S'2||)^^^, to provide information about the asym¬ 
metry of the relative velocities. Figure 18(a) indicates that the relative velocities are 
negatively skewed (Wang et al. 2000; Ray & Collins 2011). This skewness is a result of 
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Figure 18. The (a) skewness and (b) kurtosis of the relative velocities as a function of St 
for separations in the dissipation range and different values of R\. Gray filled symbols denote 
r = 1 . 7577 , and black hlled symbols denote r = 9.757], 


two contributions. First, the velocity derivatives of the underlying turbulence are neg¬ 
atively skewed, a consequence of the energy cascade (Tavoularis et al. 1978). Second, 
additional skewness arises from the path-history effect described earlier (see also Bragg 
& Collins 20146). Figure 18(a) shows by implication that at 56 ~ 1 it is the latter effect 
that dominates the skewness behavior. At even larger values of St, the effect of both 
mechanisms decreases because, with increasing Stokes number, the particle velocity dy¬ 
namics become increasingly decoupled from the small-scale fluid velocity field and their 
motion becomes increasingly ballistic in the dissipation range. 

Finally, we consider the kurtosis of the relative velocities, 'S' 4 ||/(S' 2 ||)^, in figure 18(b) 
to quantify the contributions from intermittent events in the tails of the PDFs. The 
trends are similar to those in >S'((||/(S' 2 ||)^^^, as expected, indicating that contributions 
from intermittent events become strongest for intermediate St, the smallest separations, 
and the highest Reynolds numbers. In all cases, the kurtosis is above that for a Gaussian 
distribution = 3). 


4.1.2. Inertial range relative velocity statistics 

We finally consider the inertial-range statistics of the relative velocities. In figure 12, 
we see that the relative velocities in the inertial range generally decrease with increasing 
St. This implies that the filtering mechanism (which causes the velocities to decrease 
with increasing St) dominates the path-history mechanism (which causes the velocities 
to increase with increasing St), in contrast to their relative roles in the dissipation range. 
The role reversal occurs because the path-history effect weakens as the separation is 
increased, as explained in Bragg & Collins (20146). 

For St ^ 10, the relative velocity variances appear to scale with the same scaling 
predicted by K41 for St = 0 particles. However, we observe that at St = 30, no clear 
inertial-range scaling is present. The lack of inertial scaling suggests that these particles 
are affected by their memory of large-scale turbulence throughout the entire inertial 
range. 

We now determine the scalings of the structure functions in the inertial range for 
St < 10 by computing the scaling exponents C|” and C”. Following convention (e.g., see 
Ishihara et al. 2009), we consider the scaling exponents of the relative velocity magnitudes 
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of Wy (t) and (t) here, 

W = ( li'y W oc /ii (4.8) 

and 

'5'f„iiiW = ((4-9) 

According to K41, for 77 <C r <C and St = 0, = n/3. It is well-known, how¬ 

ever, that for fluid particles, the effect of intermittency leads to a nonlinear relationship 
between and n (e.g., see Pope 2000). Kolmogorov’s refined similarity hypothesis 
(Kolmogorov 1962, hereafter ‘K62’) attempts to correct for the effect of intermittency, 
giving (for St = 0) 

C|U = ^[l-f(n-3)], (4.10) 

where 77 is typically taken to be 0.25 (Pope 2000). 

C|"j_ are shown in figure 19 at R\ = 88 and R\ = 597. For Rx = 88, we have no clear 
inertial range and therefore used extended self-similarity (Benzi et al. 1993, hereafter 
‘ESS’) to increase the scaling region for ry <C r <C i'. At Rx = 597 we have nearly 
a decade of inertial range scaling (50 r/rj < 500), and thus we can compute the 
exponents directly over this range. 

(To verify that any differences between the scaling exponents at Rx = 88 and Rx = 597 
were in fact due to Reynolds-number effects and were not merely artifacts of ESS, we 
also computed the exponents for Rx = 597 using ESS. Both methods of computing the 
exponents (directly and with ESS) gave similar results, with differences that were less 
than 8%, indicating the trends observed below are robust. We also note that while the 
inertial scaling region varies with St, we used the same fitting range for all values of St 
for consistency.) 

For St = 0, (4.10) approximates the longitudinal scaling exponents excellently for 
p < 8 at i?A = 88 (figure 19(a)), while it slightly under-predicts them at Rx = 597 
(figure 19(b)). By comparing figure 19(a) and figure 19(b), it is evident that Cy* increases 
with increasing Rx- For Rx = 88, the longitudinal scaling exponents decrease monoton- 
ically with increasing St, as was observed in Salazar & Collins (2012&). However, for 
Rx = 597, the exponents increase with St up to S't « 1 before decreasing for higher 
values of St. The reason for these trends is unclear. 

For most values of St, the transverse structure functions (figure 19(c,d)) are more 
intermittent than their longitudinal counterparts (figure 19(a,b)), in agreement with 
earlier observations (e.g., see Ishihara et al. 2009). The difference between the longitudinal 
and transverse structure functions seems to decrease as Rx increases, however, suggesting 
that it may be a low-Reynolds-number artifact (see Kerr et al. 2001; Gotoh et al. 2002; 
Shen & Warhaft 2002). 

4.2. Particle clustering 

As discussed in §1, inertial particles form clusters when placed in a turbulent flow. We 
Hrst consider a theoretical framework for understanding this clustering (§4.2.1), and then 
analyze the clustering using DNS (§4.2.2). 

4.2.1. Theoretical framework for particle clustering 

A variety of measures have been proposed to study particle clustering, including 
Voronoi' diagrams (Monchaux et al. 2010), Lyapunov exponents (Bee et al. 20066), 
Minkowski functionals (Calzavarini et al. 2008), and radial distribution functions (RDFs) 
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Figure 19. (a,b) Longitudinal and (c,d) transverse particle structure function scaling exponents 
in the inertial range for various values of St. (a,c) are for Rx — 88, and (b,d) are for Rx = 597. 
The exponents are computed from linear least-squares regression using ESS in (a,c) and directly 
in (b,d). The predicted scalings from from K41 and K62 (i.e., (4.10) with /i = 0.25) are indicated 
by the solid and dotted lines, respectively. 


(McQuarrie 1976). The RDF has distinct advantages over these other methods. The 
RDF, unlike both Minkowski functionals (Calzavarini et al. 2008) and Voronoi diagrams 
(Tagawa et al. 2012), is not biased by the number of particles simulated. Also, as Bee 
et al. (20066) noted, the accurate computation of Lyapunov exponents is numerically 
unfeasible for high-Reynolds-number simulations, while computation of the RDF is rela¬ 
tively straightforward. Finally, the RDF, unlike the other measures, has a direct relevance 
to particle collisions, since it precisely corrects the collision kernel for particle clustering 
(Sundaram & Collins 1997). 

The RDF g(r) is defined as the ratio of the number of particle pairs at a given sep¬ 
aration r to the expected number of particle pairs in a uniformly distributed particle 
field. 


9{r) 


N^/V, 

N/V ■ 


(4.11) 


Here, Ni is the number of particle pairs that lie within a shell with an average radius 
r and a radial width Ar, Vi is the volume of the shell, and N is the total number of 
particle pairs located in the total volume V. An RDF of unity corresponds to uniformly 
distributed particles, while an RDF in excess of one indicates a clustered particle field. 

Based on the findings of Bragg & Collins (2014o) we use the model of Zaichik & 
Alipchenkov (2009) as a framework for understanding the physical mechanisms govern- 
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ing particle clustering. We will validate this model against DNS data in §4.2.2. In the 
following discussion, we non-dimensionalize all variables by Kolmogorov units and use Y 
to denote the non-dimensionalized form of a variable Y. 

From Zaichik & Alipchenkov (2009), the equation describing g{f) at steady-state for 
an isotropic system is 


0 = -St (5P|+A||) 


Stg 


qP 

r02|| 


2f 


-1 


qP 

'^211 



(4.12) 


where Ay is a diffusion coefficient describing the effect of the turbulence on the dispersion 
of the particle pairs (e.g., see Bragg & Collins 2014o). We now consider (4.12) in different 
S't-regimes to consider the effect of changes in R\ within these regimes. 

In the limit St 1, (4.12) can be reduced to (see Bragg & Collins 2014a), 

0 = -r^B^iVfg - ^rg (^{5^ - {U^) , (4.13) 

where Bni is a S't-independent, non-local diffusion coefficient (see Chun et al. 2005; Bragg 
& Collins 2014o). The first term on the right-hand-side is associated with an outward 
particle diffusion which reduces clustering, while the second term on the right-hand-side 
is responsible for an inward particle drift which increases clustering. 

We therefore see that if Bni is independent of Rx, the diffusion will be independent of 
Rx- The drift is dependent on Tf^{S‘^Y — t^{TZ'^Y we see from §3.1 that ~ 

{TYy increases weakly with Rx for S't <C 1. We therefore expect the degree of clustering 
at low St to increase weakly as Rx increases. We will test this expectation against DNS 
data in §4.2.2. 

For particles with intermediate values of St, we are generally unable to simplify (4.12), 
since all terms are of comparable magnitude, and the clustering in this range is due to 
both preferential sampling and path-history effects. Bragg & Collins (2014a) showed that 
path-history effects induce an asymmetry in the particle inward and outward motions, 
causing particles to come together more rapidly than they separate, generating a net 
inward drift and increased clustering. The precise range of St over which path-history 
effects increase clustering will likely vary with Rx, but a rough guideline (based on Bragg 
& Collins 2014a) is 0.2 Y St < 0.7. Below this range, path-history effects have a negligible 
impact on particle clustering, and above this range, the path-history mechanism acts to 
diminish clustering. For the upper end of this At-range, path-history effects are the 
dominant particle-clustering mechanism (Bragg & Collins 2014a). 

We next simplify (4.12) when St > 1. As noted in §4.1.1, at sufficiently large St 
and small r/ij, the relative particle velocities are dominated by path-history effects, and 
S '211 « <S' 2 _i_. Furthermore, Ay ^ 5211 in this regime (see Bragg & Collins 2014&). Using 
these results we can simplify (4.12) in the dissipation range to the form, 

0 « -StS^^^Wfg - StgWfSP^^. (4.14) 


The overall changes in the particle clustering at high St will therefore be determined 
by the extent to which the drift coefficient (Vf-^^y) and the diffusion coefficient (^ 2 ||) 
are influenced by changes in Rx- That is, if the ratio between the drift and diffusion 
coefficients increases (decreases) with increasing Rx, the RDFs are expected to increase 
(decrease). 

We therefore take the ratio between the drift and diffusion coefficients and obtain 


V,^2"|| 




r 


(4.15) 
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r/rj r/rj 


Figure 20. RDFs for (a) low-Sf particles and (b) high-St particles at three different values of 
R\, plotted as a function of the radial separation r/rj. The Stokes numbers are indicated by the 
line labels. 

where Cy is the scaling exponent of the longitudinal relative velocity variance. (4.15) 
implies that increases (decreases) in Cy are fundamentally linked to increases (decreases) 
in the RDFs at high St. From §4.1.1, we see that Cy increases with increasing Rx for 
1 < St < 3, which suggests that g(r/rj) will increase with increasing Rx here. 

We also note that (4.14) is only applicable for high-S't particles in the dissipation 
range, and is thus unable to predict the clustering for St > 2> particles, which is primarily 
dependent on inertial-range scales. We will examine the RDFs for St > 3 from DNS data 
in §4.2.2. 

In summary, at small St, clustering may increase with increasing Rx depending upon 
whether Bni varies with Rx. Clustering at intermediate values of St will be due to both 
preferential sampling and path-history effects, though it is unclear the degree to which 
g{r/7j) will change with Rx. At high St, the degree of clustering is determined by the 
influence of path-history effects on the scaling of the relative velocity variances, which 
in turn affects the relative strengths of the drift and diffusion mechanisms. Based on our 
relative velocity data in §4.1, we expect that clustering will increase with increasing Rx 
here. We next consider DNS data to test these predictions. 

4.2.2. Particle clustering results 

In figure 20, we plot the RDFs for the different values of St considered at three different 
Reynolds numbers. Note that as the size of the simulation (and thus Rx) increases, we 
are able to calculate g{r/ri) statistics accurately at progressively smaller values of r/rj. 

In agreement with past studies (e.g., see Wang & Maxey 1993; Sundaram & Collins 
1997; Balachandar & Eaton 2010), we see that particle clustering peaks for S'! ~ 1 at all 
Reynolds numbers shown. Figure 20 also indicates that the largest particles {St > 10) 
exhibit clustering outside of the dissipation range of turbulence, and that the degree of 
clustering is independent of separation in the dissipation range. This is because large-S'l 
particles are unresponsive to the dissipative range scales and so move almost ballistically 
at these separations. The clustering that is observed for these particles is due almost 
entirely to eddies in the inertial range with timescales similar to the particle response 
time (Goto & Vassilicos 2006; Bee et al. 20106). If we make that assumption, along with 
the standard K41 approximations for the inertial range, we expect the clustering will 
depend only on e and r, and will occur at lengthscales on the order of gSt^^"^ (ElMaihy 
& Nicolleau 2005; Bee et al. 20106). We test this in figure 21 by plotting the RDFs for 
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Figure 21. RDFs for St = 20 and St = 30 particles at the three highest values of R\. The 
separations are scaled by rjSP^^ to test for inertial range scaling. The Stokes numbers are 
indicated by the line labels. 


St = 20 and St = 30 particles as a function of r/ at the three highest Reynolds 

numbers. (The two lower Reynolds numbers do not have a well-dehned inertial range, as 
noted in §2.1, and hence the above argument would not hold.) We see that the RDFs 
decrease rapidly near r/ ~ 1, suggesting that the particles are indeed clustering 

due to the influence of turbulent eddies in the inertial range with a timescale on the order 
of Tp. Refer to Bragg et al. (20156) for a recent theoretical and computational analysis 
of particle clustering in the inertial range of turbulence. 

We now discuss how the RDFs change with the Reynolds number. In §4.2.1, we argued 
that g{r/ri) might increase weakly with R\ for St 1, since —increases 

with Rx in this limit. In figure 20(a), however, we observe that g{r/ri) is essentially 
independent of Rx for St < I, which implies that the non-local correction coefficient 
Bni in (4.13) must increase weakly with Rx in a compensating way. Several authors 
have also found the level of particle clustering to be independent of Rx at small St 
(without gravity), including Collins & Keswani (2004) (from data at 65 ^ Rx ^ 152), 
Bee et al. (2007) (65 < i?A ^ 185), Bee et al. (2010a) (185 < Ra ^ 400), Ray & Collins 
(2011) (95 < Ra ^ 227), and Rosa et al. (2013) (28 < Ra ^ 304). Our data confirms 
this point up to Rx = 597. The fact that g{r/g) is independent of Rx for small Stokes 
numbers implies that the clustering mechanism is driven almost entirely by the small- 
scale turbulence, independent of any intermittency in the turbulence that occurs at higher 
Reynolds numbers. For St > 1, the RDFs increase with increasing Rx, in agreement with 
our expectations in §4.2.1. 

We note, however, that two recent studies (Onishi et al. 2013; Onishi & Vassilicos 2014) 
found that g{r/r]) decreases weakly with increasing Rx over the range 81 < Ra < 527 at 
St = 0.4 and St = 0.6. Our results do not indicate such a trend, possibly because we are 
unable to analyze g{r/g) at separations as low as those considered in Onishi et al. (2013) 
and Onishi & Vassilicos (2014). In any case, the trends with Rx at low St reported here, 
in Onishi et al. (2013) and Onishi & Vassilicos (2014), and in the rest of the literature 
are at most very weak. 

It is important to note, however, that just because g{r/g) is invariant with Rx for low- 
St particles does not necessarily imply that higher-order moments of clustering are also 
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Figure 22. Power-law fits for g{r/rj) from (4.16). (a) shows the coefficient co, and (b) shows the 
exponent ci. DNS data are shown with symbols, and the theoretical predictions from Zaichik 
& Alipchenkov (2009) (‘ZT’ and ‘ZT + DNS’), Chun et al. (2005) (‘CTl’ and ‘CT2’), and 
Gustavsson & Mehlig (2011) (‘GT’) at R\ = 597 are shown with lines and plus signs. The 
details of each of the theoretical models are discussed in the text. 


independent of R\. For example, g{r/r]) is related to the variance of the particle density 
field (Shaw et al. 2002). Higher-order moments or PDFs of the particle density field 
(e.g., see Pan et al. 2011) could also be compared at different values of R\. However, 
we found that the number of particles in our simulations was insufficient to compute 
such statistics accurately at small separations. We would likely need about an order of 
magnitude more particles to test the Reynolds-number dependence of these higher-order 
clustering moments. Refer to Yoshimoto & Goto (2007) for a more complete discussion 
on the number of particles necessary for accurate higher-order clustering statistics. 

Following Reade & Collins (2000a), we fit the RDFs by a power law of the form 

g{r/T]) « Co . (4.16) 

(Note that ci is related to the correlation dimension T >2 (Bee et al. 2007) by the relation 
Cl = 3 — T> 2 -) This allows us to compare the DNS data to several theoretical predictions 
in figure 22. For each value of R\, we computed cq and ci by fitting g{r/r]) in the range 
0.75 ^ rlrj ^ 2.75 using linear least-squares regression. For St ^ 10, we do not observe 
power-law scaling for the RDF, and thus no values of Cq and Ci are plotted here. 

To verify the arguments presented in §4.2.1, we compare the DNS values of cq and ci to 
the predicted values from Zaichik & Alipchenkov (2009) at R\ = 597. The comparisons 
are performed in two ways. In the first way (which we denote as ‘ZT’), we use the 
model of Zaichik & Alipchenkov (2009) to compute the relative velocities, and then 
use these predicted relative velocities in (4.12) to solve for the RDFs. In this manner, 
we can test the quantitative predictions of the model when no additional inputs are 
used. In the second approach (which we denote as ‘ZT -|- DNS’), we solve (4.12) with 
the particle velocities and the strain rate timescales along particle trajectories specified 
using DNS data. (The strain rate timescales are used in computing the dispersion tensor 
A. To maintain consistency in the model, we also adjusted the inertial range timescales 
through (18) in Zaichik & Alipchenkov (2003).) In both cases, we used the non-local 
diffusion correction discussed in Bragg & Collins (2014o), with Bni = 0.056. 

As expected, ‘ZT’ is only able to provide a reasonable prediction for cq and ci for 
St < 0.3. Above this point, inaccuracies in the predicted relative velocities lead to inac¬ 
curate clustering predictions, as discussed in Bragg & Collins (2014a). However, ‘ZT -|- 
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DNS’ predicts ci almost perfectly, with only slight discrepancies at St ^ 1, in agreement 
with the findings of Bragg & Collins (2014a) at a lower Reynolds number. We expect 
that these discrepancies are due to an additional drift term that was omitted in Zaichik 
& Alipchenkov (2009), as discussed in Bragg & Collins (2014a). ‘ZT + DNS’ also pro¬ 
vides reasonable predictions for cq, though the agreement is not as good as that for ci, 
possibly because cq is influenced by the inertial-range scales, which are generally more 
difficult to model. From these comparisons, we see that the model presented in §4.2.1 
is accurate, validating its use in interpreting the physical mechanisms responsible for 
particle clustering. 

We next compare our results for ci against two relations derived in Chun et al. (2005) 
in the limit of small St. The first (which we denote as ‘CTl’) uses DNS data for the strain 
and rotation rates sampled along inertial-particle trajectories to compute ci, giving 

Cl = {{sy - (ny) ■ (4.17) 

or>nl 

The second (which we denote as ‘CT2’) requires only DNS data for quantities sampled 
along fluid-particle trajectories and predicts, 
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‘CTl’ agrees well with the DNS up to St « 0.5, while ‘CT2’ only agrees well for 
St = 0.05, in agreement with Chun et al. (2005); Bragg & Collins (2014a). At higher 
values of St, both models from Chun et al. (2005) over-predict ci. As explained in Bragg 
& Collins (2014a), this over-prediction is because the theory of Chun et al. (2005) fails 
to account for the contribution of the path-history effects on the drift and diffusion 
mechanisms that govern the clustering. 

Finally, we compare our DNS values for ci against the theory from Gustavsson & 
Mehlig (2011), here denoted as ‘GT.’ The theory in Gustavsson & Mehlig (2011) predicts 
that in the limit of small r/ry, 

^^l|(xr^b (4.19) 

for n > Cl. (Note that the predictions of Zaichik & Alipchenkov (2009) and Gustavsson & 
Mehlig (2011) are equivalent when St is large, as explained in Bragg & Collins (2014a).) It 
therefore follows that for sufficiently small r/ry, ci = Cy, where is the scaling exponent 
of the relative velocity variance in the dissipation range, as computed in §4.1.1. 

We include the prediction ci = in figure 22, and see that while ‘GT’ is in excellent 
agreement with the DNS for St = 2,3, significant discrepancies exist at low St, as 
explained in Bragg & Collins (2014a). 


4.3. Collision kernel 

We now consider the kinematic collision kernel K for inertial particles, which has been 
shown to depend on both the radial distribution function and the radial relative velocities, 

K{d) = II (r = d)g(r = d), (4.20) 

where d is the particle diameter (see Sundaram & Collins 1997; Wang et al. 1998). While 
we simulate only point-particles (refer to §2.2), we compute d from St by assuming a 
given Pp/pf. To study the dependence of K(d) on Pp/pf, we consider three different 
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St St 


Figure 23. (a) The non-dimensional collision kernel K{d) as a function of St for different values 
of R\. Data are shown for pp/pf = 250 (filled black symbols), pp/pf = 1000 (open symbols), 
and pp/pf = 4000 (hlled gray symbols). Legend entries marked with f indicate data taken from 
Rosa et al. (2013) (deterministic forcing scheme, no gravity) at ppjpj = 1000. These data are 
only included in the inset, where they are compared with our results at pp/pj = 1000. (b) The 
ratio between K(d) at a given value of R\ to that at R\ = 88, to highlight any Reynolds-number 
effects for St 3. All data correspond to ppjpf = 1000. 


values for this parameter: 250, 1000, and 4000. (Note that for droplets in atmospheric 
clouds, Pp/Pf ~ 1000.) 

In general, we do not have adequate statistics to calculate g{r) or 5'!!.||(r) at r = d at 
low values of St {St < 3 for Pp/Pf = 250 and 1000, and St < 10 for PpjPf = 4000) and 
so we extrapolate from the power-law fits in §4.1.1 and §4.2.2 down to these separations, 
as was also done in Rosa et al. (2013). For larger St {St > 10 for PpjPf = 250 and 1000, 
and St > 20 for Pp/Pf = 4000), the particle diameters are sufficiently large such that 
we can compute g{d) and S'()||((i) by interpolating between data at smaller and larger 
separations. 

Following Vofikuhle et al. (2014), we compute the non-dimensional collision kernel 
K{d) = K{d)/{d^Up) = ‘iTTg{d)S^_^^{d)/up. Figure 23(a) shows K{d) for different values of 
Pp/Pf. Results from Rosa et al. (2013) (deterministic forcing scheme, no gravity, Pp/Pf = 
1000) are included in the inset to Figure 23(a). 

For St ^ 10, the collision kernels increase strongly with increasing R\^ since both the 
relative velocities and the RDFs increase with R\ here (see §4.1.1 and §4.2.2). K{d) is 
also independent of Pp/Pf here. The physical explanation is that while changes in Pp/pf 
lead to changes d, S^^^{d)/up and g{d) are largely independent of d here (see §4.1.1 and 
§4.2.2). 

Such particles, however, are generally above the size range of droplets in atmospheric 
clouds (e.g., see Ayala et al. 2008), and thus our primary focus is on the collision rates 
of smaller {St < 3) particles. K{d) is independent of Pp/Pf for 1 < S'! < 3, in agreement 
with the findings of VoBkuhle et al. (2014). In this case, while both g{d) and S^^^/up are 
dependent on d, these two quantities have opposite scalings (see §4.2.2), causing their 
product to be independent of d (and thus of pp/pf). 

For St < 3, our data show very little effect of Rx on the collision rates, and are in good 
agreement with the collision statistics from Rosa et al. (2013) at Pp/pf = 1000 (shown in 
the inset to figure 23(a)). However, since the Reynolds numbers in clouds {Rx ~ 10, 000) 
are at least an order of magnitude larger than those in the DNS, it is important to discern 
even weak trends in the collision kernel with the Reynolds number. We therefore plot the 
ratio of K{d) at a given Reynolds number to that at Rx = 88 for St ^ 3 in figure 23(b). 
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At St < 0.2, the collision statistics are almost completely independent of R\, since 
both S'^y/u^ and g are independent of R\ here (refer to §4.1.1 and §4.2.2). For larger 
St, the collision kernel very weakly decreases with increasing R\, since the mean inward 
relative velocities decrease with increasing Rx here (see §4.1.1). Finally, for 1 < At ^ 3, 
the collision kernel increases weakly as Rx increases. In this case, the increase in the RDFs 
with increasing Rx (§4.2.2) overwhelms the decrease in the relative velocities (§4.1.1), 
causing the collision kernel to increase weakly. 

These findings suggest that lower-Reynolds-number studies may in fact capture the 
essential physics responsible for droplet collisions in highly turbulent clouds. However, 
the results must be interpreted with caution for two reasons. First, the collision rates for 
St ^ S were computed by extrapolating power-law fits to very small separations, and 
it is not known if the functional form of the relative velocities and the RDFs remains 
the same at these separations. Second, even the highest Reynolds numbers in this study 
are still at least an order of magnitude smaller than those in atmospheric clouds. It is 
thus possible that the turbulence could exhibit different characteristics at much higher 
Reynolds numbers, or that the above trends in the Reynolds number, though weak, could 
lead to substantially different collision rates when Rx is increased by another order of 
magnitude. 


5. Conclusions 

We have studied the effect of particle inertia and the flow Reynolds number on particle 
dynamics at the highest Reynolds number (Rx ~ 600) and largest number of particles 
(~ 2.5 billion) to date. These simulations have provided new insights into both single- 
and two-particle statistics in homogeneous isotropic turbulence. 

We first analyzed the statistics of individual inertial particles. At large St, the particle 
motions were seen to be influenced primarily by inertial filtering. The theoretical models 
of Abrahamson (1975) and Zaichik & Alipchenkov (2008) were able to quantify the effect 
of filtering on kinetic energies and particle accelerations, respectively, in this limit, and 
provided us with a clear physical understanding of the effect of Reynolds number on 
these quantities. 

In the opposite limit (St <C 1), the particle motions were influenced primarily by 
preferential sampling, and we used the theoretical model of Chun et al. (2005) to under¬ 
stand and predict the statistics here. For St <C 1, the mean rotation rate sampled by 
the particles decreased with increasing St and Rx, since intense rotation regions became 
more prevalent and more efficient at ejecting particles (see Collins & Keswani 2004). 
As Rx increased, intense rotation regions tended to occur together with intense strain 
regions in ‘vortex sheets,’ in agreement with Yeung et al. (2012), and particles were also 
ejected from these regions, decreasing the mean strain rate sampled by the particles. In 
agreement with Salazar & Collins (20126), the particle kinetic energy increased with St 
for At <C 1 due to preferential sampling of the flow field. However, since ejections from 
vortex sheets tend to reduce the particle kinetic energy, this trend was reduced as the 
Reynolds number was increased. Fluid particle accelerations were seen to be extremely 
intermittent at high Rx, and the trends in the acceleration variance were well-captured by 
the model of Sawford et al. (2003). The particle acceleration variances decreased rapidly 
with increasing St, as inertial particles tended to be ejected from vortex tubes and vortex 
sheets, which were both characterized by very high fluid accelerations. 

We then studied the relative velocity, clustering, and collision statistics of inertial 
particles. For St 1, preferential sampling led to an increase in the longitudinal relative 
velocities and to a decrease in the transverse relative velocities, and the relative velocities 
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were generally independent of R\ for St < 0.1. At higher values of St, the particle 
motions were influenced more by path-history interactions, leading to a sharp increase in 
the relative velocities with increasing St. While the mean inward relative velocities were 
generally independent of R\ for 0.2 St < 1, the relative velocity variances increased 
weakly with increasing R\ here, a trend we attributed to either the increased scale 
separation at higher Reynolds numbers, the increased intermittency of the turbulence 
at higher Reynolds numbers, or some combination of the two. For intermediate St (1 < 
St ^ 3), the relative velocities decreased with increasing R\, which we argued was related 
to the decrease in the Lagrangian rotation timescales with increasing R\. We observed 
that the relative velocities of particles with St > 10 increased with increasing R\, since 
inertial filtering effects diminish and v!/u^ increases as the Reynolds number increases. 

We also analyzed the dissipation-range scaling exponents of the relative velocities, and 
found that particles with higher relative velocities generally had lower scaling exponents, 
since the particles were more influenced by path-history effects. Relative velocities in 
the dissipation range were seen to be strongly non-Gaussian, with the degree of non- 
Gaussianity being largest for St ~ 1, r/r; —)• 0, and high R\, suggesting that theories 
which assume a Gaussian distribution to relate the velocity variances to the mean inward 
velocities provide poor predictions for the mean inward relative velocities at particle 
contact. Higher-order inertial range structure functions were also examined and were 
observed to follow similar trends to those reported in Salazar & Collins (20126). 

We then used these trends in the relative velocities to predict the degree of clustering 
through the model of Zaichik & Alipchenkov (2009), and compared the results to DNS 
data. The trends in the RDFs at low St were tied to preferential sampling effects, which 
increased the inward particle drift, as was found in Chun et al. (2005). The RDFs were 
independent of R\ here, in agreement with Collins & Keswani (2004); Ray & Collins 
(2011); Rosa et al. (2013), suggesting that the non-local coefhcient Bni (see Chun et al. 
2005; Bragg & Collins 2014a) must weakly increase as R\ increases. (We were unable to 
test higher-order measures of clustering to determine if they were affected by changes in 
R\ due to the limitations in the number of particles that could be simulated.) 

At high St, the degree of clustering was tied to the influence of path-history effects on 
the particle drift and diffusion, as explained in Bragg & Collins (2014a). By simplifying 
the model of Zaichik & Alipchenkov (2009) in this limit, we showed that changes in the 
scaling exponents of the relative velocity variances directly affected the drift and diffu¬ 
sion mechanisms, which in turn altered the clustering. The scaling exponents generally 
increased with increasing R\ (suggesting that path-history effects became less impor¬ 
tant), which in turn led to increased levels of clustering. For St ^ 10 and R\ ^ 224, 
particles were seen to cluster in the inertial range of turbulence, and the separation at 
which clustering decreased was predicted accurately by inertial-range scaling arguments. 

For St < 3, the RDFs exhibited power-law scaling, consistent with Reade & Collins 
(2000a). The full model of Zaichik & Alipchenkov (2009) (without any inputs from 
the DNS) was able to predict the power-law coefficient cq and power-law exponent ci 
accurately only for St < 0.4 due to errors in the predicted relative velocities. However, 
when these relative velocities (and the associated Lagrangian timescales) were specified 
from the DNS, the model in Zaichik & Alipchenkov (2009) provided excellent predictions 
for Cl and reasonable predictions for cq, as was also found in Bragg & Collins (2014a) at 
a lower Reynolds number. We also tested the DNS against two model predictions from 
Chun et al. (2005), one which required only fluid particle statistics from the DNS, and 
one which required strain and rotation statistics along particle trajectories. The former 
prediction was in acceptable agreement with the DNS only for St = 0.05, while the latter 
prediction was in good agreement up to St « 0.5, in agreement with Chun et al. (2005); 
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Bragg & Collins (2014o). Finally, we found that the theory of Gustavsson & Mehlig 
(2011) was able to predict Ci well for St = 2 and St = 3. 

We used the relative velocity and RDF data to compute the kinematic collision kernel 
for inertial particles (Sundaram & Collins 1997), and found that this quantity varied 
only slightly with Reynolds number (under 50% when Rx changed by a factor of 7) for 
0 ^ St ^ S. Our collision kernels were in good agreement with those computed by Rosa 
et al. (2013). 

As mentioned in §1, one of the primary motivations for this study was to determine the 
extent to which turbulence-induced collisions are responsible to the rapid growth rate of 
droplets observed in warm, cumulus clouds. Our observations indicate that the collision 
rates of like particles are generally unaffected by changes in the Reynolds number, which 
suggests that relatively low-Reynolds-number simulations may allow us to study the es¬ 
sential physics of droplet collisions in highly turbulent atmospheric clouds. One promising 
avenue of future work would be to determine the droplet growth rates predicted by these 
collision kernels, either by solving an associated kinetic equation (Xue et al. 2008; Wang 
& Grabowski 2009) or by simulating the particle collision and coalescence process directly 
(Reade & Collins 20006). 

Finally, we note that it is unclear to what extent these conclusions would be altered 
if gravity were incorporated in the particle dynamics, since the introduction of gravity 
will likely cause particles to preferentially sample certain regions of the flow, and will 
alter the residence time of particles around certain flow features (e.g., see Wang & Maxey 
1993; Davila & Hunt 2001; Good et al. 2014). We will analyze the effect of gravity on 
inertial particle motion in turbulence in Part II (Ireland et al. 2015). 


Acknowledgements 

The authors gratefully acknowledge Parvez Sukheswalla for helpful discussions regard¬ 
ing this work. This work was supported by the National Science Foundation through 
CBET grants 0756510 and 0967349, and through a graduate research fellowship awarded 
to PJI. Additional funding was provided by Cornell University. We would also like to ac¬ 
knowledge high-performance computing support from Yellowstone (ark:/85065/d7wd3xhc) 
provided by NCAR’s Computational and Information Systems Laboratory through grants 
ACOROOOl and P3509I057, sponsored by the National Science Foundation. 


REFERENCES 

Abrahamson, J. 1975 Collision rates of small particles in a vigorously turbulent fluid. Chem. 
Eng. Sci. 30 , 1371-1379. 

Ashurst, W. T., Kerstein, A. R., Kerr, R. M. & Gibson, C. H. 1987 Alignment ofvorticity 
and scalar gradient with strain rate in simulated Navier-Stokes turbulence. Phys. Fluids 
30 (8), 2343-2353. 

Ayala, O., Rosa, B., Wang, L.-P. & Grabowski, W. W. 2008 Effects of turbulence on 
the geometric collision rate of sedimenting droplets, part 1. results from direct numerical 
simulation. New J. Phys. 10 , 075015. 

Ayyalasomayajula, S., Warhaft, Z. & Collins, L. R. 2008 Modeling inertial particle ac¬ 
celeration statistics in isotropic turbulence. Phys. Fluids 20 , 094104. 

Balachandar, S. & Eaton, J. K. 2010 Turbulent dispersed multiphase flow. Annu. Rev. Fluid 
Meeh. 42 , 111-133. 

Beg, j., Biferale, L., Boffetta, G., Celani, A., Cencini, M., Lanotte, A. S., Musac- 
CHIO, S. & Toschi, F. 2006o Acceleration statistics of heavy particles in turbulence. J. 
Fluid Mech. 550 , 349-358. 



34 


P. J. Ireland et al. 


Bec, J., Biferale, L., Boffetta, G., Cencini, M., Musacchio, S. & Toschi, F. 2006b 
Lyapunov exponents of heavy particles in turbulence. Phys. Fluids 18 , 091702. 

Bec, j., Biferale, L., Cencini, M., Lanotte, A. S., Musacchio, S. & Toschi, F. 2007 
Heavy particle concentration in turbulence at dissipative and inertial scales. Phys. Rev. 
Lett. 98 , 084502. 

Bec, j., Biferale, L., Cencini, M., Lanotte, A. S. & Toschi, F. 2010a Intermittency in 
the velocity distribution of heavy particles in turbulence. J. Fluid Mech. 646 , 527-536. 

Bec, j., Biferale, L., Lanotte, A. S., Scagliarini, A. & Toschi, F. 2010b Turbulent pair 
dispersion of inertial particles. J. Fluid Mech. 645 , 497-528. 

Benzi, R., Ciliberto, S., Tripiccione, R., Baudet, C., Massaioli, F. & Succi, S. 1993 
Extended self-similarity in turbulent flows. Phys. Rev. E 48 , R29-R32. 

Biferale, L., Boffetta, G., Celani, A., Lanotte, A. & Toschi, F. 2005 Particle trapping 
in three-dimensional fully developed turbulence. Phys. Fluids 17 , 021701. 

Bragg, A. D. & Collins, L. R. 2014a New insights from comparing statistical theories for 
inertial particles in turbulence: 1. Spatial distribution of particles. New J. Phys. 16 , 055013. 

Bragg, A. D. & Collins, L. R. 2014b New insights from comparing statistical theories for 
inertial particles in turbulence: II: Relative velocities. New J. Phys. 16 , 055014. 

Bragg, A. D., Ireland, P. J. & Collins, L. R. 2015a Forward and backward in time 
dispersion of fluid and inertial particles in isotropic turbulence. Phys. Fluids Submitted, 
eprint arXiv:1403.5502. 

Bragg, A. D., Ireland, P. J. & Collins, L. R. 2015b Mechanisms for the clustering of 
inertial particles in the inertial range of isotropic turbulence. Phys. Rev. E In review, 
eprint arXiv:1411.7422. 

Calzavarini, E., Kersgher, M., Lohse, D. & Toschi, F. 2008 Dimensionality and morphol¬ 
ogy of particle and bubble clusters in turbulent flow. J. Fluid Meeh. 607 , 13-24. 

Chun, J., Koch, D. L., Rani, S., Ahluwalia, A. & Collins, L. R. 2005 Clustering of aerosol 
particles in isotropic turbulence. J. Fluid Mech. 536 , 219-251. 

Collins, L. R. & Keswani, A. 2004 Reynolds number scaling of particle clustering in turbulent 
aerosols. New J. Phys. 6, 119. 

Computational and Information Systems Laboratory 2012 Yellowstone: IBM iDataPlex 
System (University Community Computing). http://n2t.net/ark:/85065/d7wd3xhc. 

Cuzzi, J. N., Hogan, R. C., Paque, J. M. & Dobrovolskis, A. R. 2001 Size-selective 
concentration of chrondrules and other small particles in protoplanetary nebula turbulence. 
Astrophysieal J. 546 , 496-508. 

Davila, J. & Hunt, J. C. R. 2001 Settling of small particles near vortices and in turbulence. 
J. Fluid Mech 440 , 117-145. 

Devenish, B. j., Bartello, P., Brenguier, J.-L., Collins, L. R., Grabowski, W. W., 
IJzermans, R. H. a., Malinowski, S. P., Reeks, M. W., Vassilicos, J. C., Wang, 
L.-P. & Warhaft, Z. 2012 Droplet growth in warm turbulent clouds. Q. J. R. Meteorol. 
Soc. 138 , 1401-1429. 

Durham, W. M., Climent, E., Barry, M., Lillo, E. D., Boffetta, G., Cencini, M. & 
Stocker, R. 2013 Turbulence drives microscale patches of motile phytoplankton. Nat. 
Commun. 4 (2148), 1-7. 

Eaton, J. K. & Eessler, J. R. 1994 Preferential concentration of particles by turbulence. Int. 
J. Multiphase Flow 20 , 169-209. 

Elghobashi, S. E. & Truesdell, G. C. 1992 Direct simulation of particle dispersion in a 
decaying isotropic turbulence. J. Fluid Mech. 242 , 655. 

Elghobashi, S. E. & Truesdell, G. C. 1993 On the two-way interaction between homoge¬ 
neous turbulence and dispersed particles, i: Turbulence modification. Phys. Fluids A 5 , 
1790-1801. 

ElMaihy, a. & Nicolleau, E. 2005 Investigation of the dispersion of heavy-particle pairs and 
Richardson’s law using kinematic simulation. Phys. Rev. E 71 , 046307. 

Falkovich, G., Fouxon, A. & Stepanov, M. G. 2002 Acceleration of rain initiation by cloud 
turbulence. Nature 419 , 151-154. 

Falkovich, G. & Pumir, A. 2007 Sling effect in collisions of water droplets in turbulent clouds. 
J. Atm. Sci. 64 , 4497. 

Good, G. H., Ireland, P. J., Bewley, G. P., Bodenschatz, E., Collins, L. R. & 



Reynolds-numher effects on inertial particle dynamics, Part I. 


35 


Warhaft, Z. 2014 Settling regimes of inertial particles in isotropic turbulence. J. Fluid 
Mech. 759, R3. 

Goto, S. & Vassilicos, J. C. 2006 Self-similar clustering of inertial particles and zero- 
acceleration points in fully developed two-dimensional turbulence. Phys. Fluids 18, 115103. 

Gotoh, T., Fukayama, D. & Nakano, T. 2002 Velocity field statistics in homogeneous steady 
turbulence obtained using a high-resolution direct numerical simulation. Phys. Fluids 14, 
1065-1081. 

Grabowski, W. W. & Wang, L.-P. 2013 Growth of cloud droplets in a turbulent environment. 
Annu. Rev. Fluid Mech. 45, 293-324. 

Gustavsson, K. & Mehlig, B. 2011 Distribution of relative velocities in turbulent aerosols. 
Phys. Rev. E 84, 045304. 

Hill, R. J. 2002 Scaling of acceleration in locally isotropic turbulence. J. Fluid Mech. 452, 
361-370. 

VAN Hinsberg, M. a. T., ten True Bookkkamp, J. H. M., Toschi, F. & Glergx, H. 
J. H. 2013 Optimal interpolation schemes for particle tracking in turbulence. Phys. Rev. 
E 87, 043307. 

IJzermans, R. H. a., Meneguz, E. & Reeks, M. W. 2010 Segregation of particles in in¬ 
compressible random flows: singularities, intermittency and random uncorrelated motion. 
J. Fluid Mech. 653, 99-136. 

Ireland, P. J., Bragg, A. D. & Gollins, L. R. 2015 The effect of Reynolds number on 
inertial particle dynamics in isotropic turbulence. Part II: Simulations with gravitational 
effects. J. Fluid Mech. Submitted. 

Ireland, P. J., Vaithianathan, T., Sukheswalla, P. S., Ray, B. & Gollins, L. R. 2013 
Highly parallel particle-laden flow solver for turbulence research. Comput. Fluids 76, 170- 
177. 

ISHIHARA, T., Gotoh, T. & Kaneda, Y. 2009 Study of high-Reynolds-number isotropic tur¬ 
bulence by direct numerical simulation. Annu. Rev. Fluid Mech. 41, 165-180. 

ISHiHARA, T., Kaneda, Y., Yokokawa, M., Itakura, K. & Uno, A. 2007 Small-scale statis¬ 
tics in high-resolution direct numerical simualtion of turbulence: Reynolds number depen¬ 
dence of one-point velocity gradient statistics. J. Fluid Mech. 592, 335-366. 

Kaneda, Y., Ishihara, T., Yokokawa, M., Itakura, K. & Uno, A. 2003 Energy dissipation 
rate and energy spectrum in high resolution direct numerical simulations of turbulence in 
a periodic box. Phys. Fluids 15, L21-L24. 

Kerr, R. M., Meneguzzi, M. & Gotoh, T. 2001 An inertial range crossover in structure 
functions. Phys. Fluids 13, 1985-1994. 

Kolmogorov, A. N. 1941 The local structure of turbulence in an incompressible viscous fluid 
for very large Reynolds numbers. Dokl. Akad. Nauk. SSSR 30, 299-303. 

Kolmogorov, A. N. 1962 A refinement of previous hypotheses concerning the local structure 
of turbulence in a viscous incompressible fluid at high Reynolds number. J. Fluid Mech. 
13, 82-85. 

Maxey, M. R. 1987 The motion of small spherical particles in a celluar flow field. Phys. Fluids 
30, 1915-1928. 

Maxey, M. R. & Riley, J. J. 1983 Equation of motion for a small rigid sphere in a nonuniform 
flow. Phys. Fluids 26, 883-889. 

MgQuarrie, D. a. 1976 Statistical Mechanics. Harper & Row, New York. 

Meneveau, C. 2011 Lagrangian dynamics and models of the velocity gradient tensor in turbu¬ 
lent flows. Annu. Rev. Fluid Mech. 43, 219-245. 

Monchaux, R., Bourgoin, M. & Cartellier, A. 2010 Preferential concentration of heavy 
particles: A Voronoi analysis. Phys. Fluids 22, 103304. 

Onishi, R., Takahashi, K. & Vassilicos, J. G. 2013 An efficient parallel simulation of in¬ 
teracting inertial particles in homogeneous isotropic turbulence. J. Comput. Phys. 242, 
809-827. 

Onishi, R. & Vassilicos, J. G. 2014 Gollision statistics of inertial particles in two-dimensional 
homogeneous isotropic turbulence with an inverse cascade. J. Fluid Mech. 745, 279-299. 

Orszag, S. a. & Patterson, G. S. 1972o Numerical simulation of three-dimensional homo¬ 
geneous isotropic turbulence. Phys. Rev. Lett. 28, 76-79. 



36 


P. J. Ireland et al. 


Orszag, S. a. & Patterson, G. S. 19726 Numerical simulation of turbulence. Springer-Verlag, 
New York. 

Pan, L. & Padoan, P. 2010 Relative velocity of inertial particles in turbulent flows. J. Fluid 
Meeh. 661, 73-107. 

Pan, L. & Padoan, P. 2013 Turbulence-indnced relative velocity of dnst particles i: identical 
particles. ApJ 776, 12. 

Pan, L., Padoan, P., Scalo, J., Kritsuk, A. G. & Norman, M. L. 2011 Turbulent clustering 
of protoplanetary dust and planetesimal formation. ApJ 740, 6. 

Pekurovsky, D. 2012 P3DFFT: A framework for parallel computations of Fourier transforms 
in three dimensions. SIAM J. Sci. Comput. 34 (4), C192-C209. 

Pope, S. B. 2000 Turbulent Flows. Cambridge University Press, New York. 

Pruppacher, H. R. & Klett, J. D. 1997 Microphysics of Clouds and Precipitation. Kluwer, 
Dordrecht. 

Ray, B. & Collins, L. R. 2011 Preferential concentration and relative velocity statistics of 
inertial particles in Navier-Stokes turbulence with and without hltering. J. Fluid Mech. 
680, 488-510. 

Ray, B. & Collins, L. R. 2013 Investigation of sub-kolmogorov inertial particle pair dynamics 
in turbulence using novel satellite particle simulations. J. Fluid Mech. 720, 192-211. 

Reade, W. C. & Collins, L. R. 2000a Effect of preferential concentration on turbulent collision 
rates. Phys. Fluids 12, 2530-2540. 

Reade, W. C. & Collins, L. R. 20006 A numerical study of the particle size distribution of 
an aerosol undergoing turbulent coagulation. J. Fluid Mech. 415, 45-64. 

Rosa, B., Parishani, H., Ayala, O., Grabowski, W. W. & Wang, L. P. 2013 Kinematic 
and dynamic collision statistics of cloud droplets from high-resolution simulations. New J. 
Phys. 15, 045032. 

Salazar, J. P. L. C. & Collins, L. R. 2012o Inertial particle acceleration statistics in tur¬ 
bulence: effects of filtering, biased sampling, and flow topology. Phys. Fluids 24, 083302. 

Salazar, J. P. L. C. & Collins, L. R. 20126 Inertial particle relative velocity statistics in 
homogeneous isotropic turbulence. J. Fluid Mech. 696, 45-66. 

Sawford, B. L., Yeung, P.-K., Borgas, M. S., La Porta, P. V. A., Crawford, A. M. & 
Bodenschatz, E. 2003 Conditional and unconditional acceleration statistics in turbulence. 
Phys. Fluids 15, 3478-3489. 

Shaw, R. A. 2003 Particle-turbulence interactions in atmospheric clouds. Annu. Rev. Fluid 
Mech. 35, 183-227. 

Shaw, R. A., Kostinski, B. & Larsen, M. L. 2002 Towards quantifying droplet clustering in 
clouds. Q. J. R. Meteorol. Soc. 128, 1043-1057. 

Shen, X. & Warhaft, Z. 2002 Longitudinal and transverse structure functions in sheared and 
unsheared wind-tunnel turbulence. Phys. Fluids 14, 370-381. 

SiEBERT, H., Lehmann, K. & Wendisgh, M. 2006 Observations of small-scale turbulence and 
energy dissipation rates in the cloudy boundary layer. J. Atmos. Sci. 63, 1451-1466. 

Soria, J., Sondergaard, R., Cantwell, B. J., Chong, M. S. & Perry, A. E. 1994 A study 
of the fine-scale motions of incompressible time-developing mixing layers. Phys. Fluids 6 (2), 
871-884. 

Spalart, P. R. 1988 Direct simulation of a turbulent boundary layer up to Rg — 1410. J. Fluid 
Mech. 187, 61-98. 

Squires, K. D. & Eaton, J. K. 1991 Preferential concentration of particles by turbulence. 
Phys. Fluids A 3, 1169-1178. 

Sundaram, S. & Collins, L. R. 1997 Collision statistics in an isotropic, particle-laden turbu¬ 
lent suspension 1. Direct numerical simulations. J. Fluid Mech. 335, 75-109. 

Sundaram, S. & Collins, L. R. 1999 A numerical study of the modulation of isotropic tur¬ 
bulence by suspended particles. J. Fluid Mech. 379, 105-143. 

Tagawa, Y., Mergado, j. M., Prakash, V. N., Calzavarini, E., Sun, C. & Lohse, D. 
2012 Three-dimensional Lagrangian Voronoi analysis for clustering of particles and bubbles 
in turbulence. J. Fluid Mech. 693, 201-215. 

Tavoularis, S., Bennett, J. C. & Corrsin, S. 1978 Velocity-derivative skewness in small 
Reynolds number, nearly isotropic turbulence. J. Fluid Mech. 88, 63-69. 

VAN Hinsberg, M. a. T., Thije Boonkkamp, j. H. M., Tosghi, F. & Clercx, H. j. H. 



Reynolds-number effects on inertial particle dynamics, Part I. 


37 


2012 On the efficiency and accuracy of interpolation methods for spectral codes. SIAM J. 
Sci. Comput. 34 (4), B479-B498. 

VOSSKUHLE, M., PuMiR, A., Leveque, E. & WiLKiNSON, M. 2014 Prevalence of the sling 
effect for enhancing collision rates in turbulent suspensions. J. Fluid Mech. 749, 841-852. 

VOTH, G. A., La Porta, A., Crawford, A. M., Alexander, J. & Bodenschatz, E. 2002 
Measurement of particle accelerations in fully developed turbulence. J. Fluid Mech. 469, 
121-160. 

Wang, L.-P. & Grabowski, W. W. 2009 The role of air turbulence in warm rain initiation. 
Atmos. Sci. Let. 10, 1-8. 

Wang, L.-P. & Maxey, M. R. 1993 Settling velocity and concentration distribution of heavy 
particles in homogeneous isotropic turbulence. J. Fluid Mech. 256, 27-68. 

Wang, L.-P., Wexler, A. S. & Zhou, Y. 1998 Statistical mechanical descriptions of turbulent 
coagulation. Phys. Fluids 10, 2647-2651. 

Wang, L.-P., Wexler, A. S. & Zhou, Y. 2000 Statistical mechanical description and modeling 
of turbulent collision of inertial particles. J. Fluid Mech. 415, 117-153. 

Wilkinson, M. & Mehlig, B. 2005 Caustics in turbulent aerosols. Europhys. Lett. 71, 186-192. 

Wilkinson, M., Mehlig, B. & Bezuglyy, V. 2006 Caustic activation of rain showers. Phys. 
Rev. Lett. 97, 048501. 

WiTKOWSKA, A., Brasseur, J. G. & JuvE, D. 1997 Numerical study of noise from isotropic 
turbulence. J. Comput. Acoust. 5, 317-336. 

XuE, Y., Wang, L.-P. & Grabowski, W. W. 2008 Growth of cloud droplets by turbulent 
collision-coalescence. J. Atmos. Sci. 65, 331-356. 

Yeung, P. K., Donzis, D. A. & Sreenivasan, K. R. 2012 Dissipation, enstrophy, and pressure 
statistics in turbulence simulations at high reynolds numbers. J. Fluid Mech. 700, 5-15. 

Yeung, P. K. & Pope, S. B. 1989 Lagrangian statistics from direct numerical simulations of 
isotropic turbulence. J. Fluid Mech. 207, 531-586. 

Yeung, P. K., Pope, S. B., Lamorgese, A. G. & Donzis, D. A. 2006 Acceleration and 
dissipation statistics of numerically simulated isotropic turbulence. Phys. Fluids 18 (6), 
065103. 

Yoshimoto, H. & Goto, S. 2007 Self-similar clustering of inertial particles in homogeneous 
turbulence. J. Fluid Mech. 577, 275-286. 

Yudine, M. I. 1959 Physical considerations on heavy-particle dispersion. Adv. Geophys. 6, 
185-191. 

Zaichik, L. I. & Alipchenkov, V. M. 2003 Pair dispersion and preferential concentration of 
particles in isotropic turbulence. Phys. Fluids 15, 1776-1787. 

Zaichik, L. I. & Alipchenkov, V. M. 2008 Acceleration of heavy particles in isotropic tur¬ 
bulence. Int. J. Multiphase Flow 34 (9), 865-868. 

Zaichik, L. I. & Alipchenkov, V. M. 2009 Statistical models for predicting pair dispersion 
and particle clustering in isotropic turbulence and their applications. New J. Phys. 11, 
103018. 

Zaichik, L. L, Simonin, O. & Alipchenkov, V. M. 2003 Two statistical models for predicting 
collision rates of inertial particles in homogeneous isotropic turbulence. Phys. Fluids 15, 
2995-3005. 



